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ABSTRACT 


Studying  seismic  wave  propagation  by  ray  techniques  has 
been  popular  in  the  geophysical  community  for  some  time.  Due 
to  the  low  cost  of  computation  and  the  ability  to  associate 
individual  events  on  synthetic  seismograms  with  particular 
paths  of  energy  propagation  (ray  paths),  asymptotic  ray 
theory  is  one  of  the  most  practical  of  the  ray  methods. 
Although  the  theoretical  aspects  of  this  theory  have  been 
well  documented,  much  greater  difficulties  of  a 
computational  nature  must  be  overcome  to  truly  establish 
this  theory.  This  thesis  addresses  itself  to  three  of  these 
topics.  First,  the  role  of  the  non-leading  term  in  the  ray 
series  in  inhomogeneous  media  is  examined.  Second,  the 
accuracy  of  the  ray  method  for  reflected  and  head  waves  is 
checked  by  comparison  with  direct  numerical  integration  of 
the  displacement  potentials.  Finally,  several  chapters  are 
devoted  to  the  development  and  testing  of  computational 
methods  applicable  to  numerical  modelling  of  seismic  body 
waves  in  laterally  inhomogeneous  media  with  curved 
interfaces.  Travel  times,  amplitude-distance  curves  and 
synthetic  seismograms  are  computed  for  several  models. 
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1 .  GENERAL  INTRODUCTION 


Elastic  wave  propagation  problems  are  of  paramount 
importance  in  all  branches  of  seismology.  For  complex  media, 
analytical  solutions  to  the  equations  of  motion  are  unknown 
and  seismologists  must  turn  to  numerical  techniques  to 
investigate  wave  fields.  One  common  approach  is  a  direct 
computational  solution  to  the  elastodynamic  equation.  This 
technique  generally  has  some  variation  of  the  finite 
difference  algorithm  at  its  hub.  If  properly  applied,  this 
method  will  give  a  total  solution  including  refracted  and 
reflected  waves,  surface  waves,  and  diffracted  waves. 
Unfortunately,  on  the  resulting  synthetic  seismograms  no 
identification  of  particular  phases  can  be  made.  For  this 
purpose,  one  can  use  ray  methods.  This  thesis  will  be 
devoted  to  various  computational  topics  related  to  one  of 
these  methods,  namely  asymptotic  ray  theory,  often  termed 
the  ray  series  method  in  European  literature. 


The  basic  tenets  of  asymptotic  ray  theory  were 
developed  simultaneously  and  independently  by  Babich  and 
Alekseev  (1958)  in  the  Soviet  Union  and  Karal  and  Keller 
(1959)  in  America.  This  ray  theory  is  powerful  in  that  it 
can  be  applied  to  the  elastodynamic  equation  in  general 
inhomogeneous  media.  However,  in  order  to  fully  employ 
asymptotic  ray  theory,  some  very  lengthy  computations 
related  to  the  seismic  body  wave  field  must  be  performed. 
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This  thesis  presents  techniques  and  numerical  results 
for  three  computational  topics  related  to  asymptotic  ray 
theory.  The  basic  principles  of  the  theory  are  presented  in 
Chapter  2.  For  a  particular  choice  of  inhomogeneous  media, 
the  second  term  in  the  ray  series  (which  is  inversely 
proportional  to  frequency)  is  evaluated  numerically  and 
compared  with  the  leading  term  which  coincides  with 
geometrical  ray  theory.  The  most  interesting  feature  of  the 
non- leading  term  is  that  for  P  waves,  it  has  components  both 
parallel  to  and  perpendicular  to  the  ray.  Due  to  the 
frequency  dependence  of  this  term,  it  will  be  shown  that  it 
is  of  little  significance  relative  to  the  leading  term  for 
source-receiver  distances  less  than  the  critical  distance 
and  for  small  degrees  of  inhomogeneity  in  the  earth. 


In  Chapter  3,  the  leading  term  of  the  ray  series  is 
compared  with  direct  numerical  integration  of  the  expression 
representing  displacement.  In  the  neighborhood  of  a  critical 
point,  where  ray  theory  is  not  valid,  the  computations  are 
compared  with  a  more  exact  wave  method.  Some  results  of  this 
work  have  been  published  by  Marks  and  Hron  (1978). 


The  last  topic,  numerical  modelling  of  seismic  body 
waves  in  laterally  inhomogeneous  media  with  curved 
interfaces,  is  discussed  in  Chapters  4-7.  Chapter  4  deals 
with  methods  of  ray  tracing  which  are  currently  in 
existence.  These  techniques  involve  either  the  solution  of  a 
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set  of  differential  equations  at  many  points  along  the  ray 
or  the  inversion  of  a  large  matrix.  Due  to  the  lack  of 
sufficient  computer  time  to  properly  utilize  these 
approaches,  an  approximation  is  presented  in  Chapter  5  and 
compared  with  the  other  techniques.  Termed  the  circular 
approximation,  it  concerns  itself  only  with  the 
inhomogeneities  of  the  medium  in  the  vicinity  of  the 
raypath.  With  a  transformation  of  coordinate  systems,  one  is 
allowed  to  write  compact  analytical  expressions  for  the  ray 
paths.  This  produces  a  tremendous  saving  in  computing  time 
and  allows  one  to  calculate  economically,  yet  accurately 
values  of  travel  times,  amplitudes  and  finally  synthetic 
seismograms . 


The  works  in  Chapters  4-5  are  in  essence,  a  paper  which 
was  read  to  the  European  Sei smo logical  Commission  (Marks  and 
Hron  1978a) . 


One  prerequisite  for  computing  ray  amplitude  is  to 
calculate  the  geometrical  spreading  of  the  wave  front. 
Analytic  techniques  which  perform  this  task  are  generally 
cast  as  sets  of  differential  equations  or  time  integration 
approaches.  These,  however,  consume  vast  amounts  of  computer 
time  and  are  better  replaced  by  a  modification  of  the 
circular  approximation  technique.  Chapter  6  describes  this 
alteration  as  well  as  two  of  the  analytic  means  of 


determining  wave  front  spreading.  A  number  of  computational 


. 
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examples  and  comparisons  are  presented  which  indicate  the 
expediency  of  this  approximation. 

Chapter  7  concerns  itself  with  a  general  procedure 
linking  together  the  seismic  ray  tracing  techniques  and 
geometrical  spreading  of  wave  fronts  to  produce  synthetic 
seismograms  for  laterally  varying  media.  The  circular 
approximation  method  is  the  basis  for  most  work  performed  in 
Chapter  7.  The  work  presented  in  Chapter  7  is  the  basis  for 
a  paper  presented  to  the  Society  of  Exploration 
Geophysicists  (Marks  and  Hron  1979b). 


2.  ACCURACY  OF  GEOMETRICAL  RAY  THEORY  IN  SEISMOLOGY 


2 . 1  Introduction 


Applied  seismology  is  concerned  with  locating  and 
geologically  interpreting  structures  within  the  earth.  Rays 
are  often  used  as  a  means  of  following  the  path  that  seismic 
wave  energy  takes  in  propagating  from  source  to  receiver.  In 
short,  the  laws  of  geometrical  ray  theory  provide  the 
groundwork  for  many  of  the  existing  seismic  i nterpretat ion 
methods.  It  is  rather  surprising  to  notice  how  uncritically 
the  approximation  is  accepted.  In  this  chapter  we  shall  show 
results  of  a  test  of  the  accuracy  of  geometrical  ray  theory. 

What  is  meant  by  geometrical  ray  theory  in  seismology? 
Basically,  it  is  the  concept  that  seismic  energy  travels 
through  an  elastic  medium  along  well-defined  paths  termed 
rays.  The  geometry  of  these  paths  is  governed  by  the  well 
known  Snell's  law.  This  definition,  however,  gives  no 
information  on  energy  properties  (amplitudes).  For  this,  it 
is  necessary  to  turn  to  the  law  of  conservation  of  energy. 
This  states  that  the  energy  flux  across  the  wave  front  must 
be  equal  at  all  times  as  long  as  the  elastic  parameters  of 
the  medium  are  continuous.  Thus,  if  a  value  for  the 
displacement  is  known  at  some  reference  point  on  the  ray,  it 
is  a  simple  matter  to  compute  displacement  at  any  other 
point  along  the  ray  by  comparing  surface  areas  of  the  wave 
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fronts,  velocities  and  densities  at  these  two  points  plus 
the  effect  of  energy  partitioning  at  interfaces.  These 
concepts  will  be  discussed  in  more  detail  later  on.  For  the 
remainder  of  this  chapter  the  term  "geometrical  ray  theory" 
will  imply  both  the  energy  conservation  law  and  Snell's  Law. 

Asymptotic  ray  theory  will  be  applied  to  the 
propagation  of  elastic  waves  through  an  inhomogeneous, 
perfectly  elastic  and  isotropic  medium.  The  density  of  the 


medium  is  denoted  p  and  its  Lame  parameters  m  and 

A  .  The  el astodynami c  equation  for  particle  displacement, 

— ^ 

w  ,  i s  given  by 


(2.1) 


+  Vyx(VxW)  +  2  ( VM'  •  V )  W 


where  t  is  time.  Using  the  basic  tenets  as  set  forth  by  Hron 
and  Kanasewich  (1971)  or  Cerveny  and  Ravindra  (1971)  a  time 
harmonic  solution  to  2.1  is  assumed  to  be  given  by  a  series 
in  inverse  powers  of  frequency,  w 


W  =  exp  (io)  (t-x  (x,y,  z)  )  T 

k=0 


CO 


(2.2) 


This  expansion  is  termed  a  ray  series  with  t 
function  and 


the  phase 


_ 
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K  ,  (k=0 , 1 , 2 , . . . )  the  amplitude  coefficients  of  the 

ser i es . 


In  this  chapter  we  attempt  quantitively  to  describe  the 
frequency  dependent  connections  between  the  magnitudes  of 
the  first  two  terms  in  2.2  and  the  degree  of  inhomogeneity 
of  the  medium.  These  results  will  also  illustrate  the  error 
which  arises  when  only  the  leading  term 

w  =  exp  (iui  (t-x)  )i$Q  (2.3) 


is  taken.  This  expression,  termed  the  zeroth  order 
approximation,  is  that  which  can  be  derived  through  the  laws 
of  geometrical  ray  theory.  On  the  other  hand  we  term 


w 


exp  ( iu)  ( t-T  )  )  (W  +S  /iw) 


(2.4) 


to  be  the  first  order  approximation  which  essentially  serves 


as  a  correction  represented  by 


- y 

w. 


elu(t-T)/iw 


Let  us  briefly  mention  two  properties  of  the  first 
order  approximation  as  opposed  to  the  zeroth  order.  The 
first  of  these  is  the  loss  of  polarization  of  particle 
displacement.  W1  ,  in  general  has  components  normal  to  the 
wave  front,  corresponding  to  compressiona 1  wave  motion,  and 
tangent  to  the  wave  front,  corresponding  to  shear  wave 
motion.  For  P  waves,  the  component  normal  to  the  wave  front 
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is  termed  the  principal  component  while  the  component 
tangent  to  the  wave  front  is  called  the  additional 
component.  Analogous  definitions  exist  for  S  waves. 
Secondly,  the  correction  term  in  2.4  is  frequency  dependent 
implying  justification  of  the  zeroth  order  approximation  for 
high  frequency  sources. 

A  perusal  of  seismological  literature  demonstrates  that 
little  work  has  been  performed  on  determining  the  magnitude 
of  any  term  in  the  ray  series  other  than  the  zeroth  order 
approximation.  Brekhovskikh  (1960)  carried  out  an  analysis 
of  the  field  of  a  reflected  spherical  electromagnetic  wave 
as  well  as  an  acoustic  wave.  However,  the  approach  used  in 
that  book  is  inappropriate  for  solving  the  more  intricate 
el astodynamic  equation.  Pod'  '  yapol '  skiy  (1966)  calculated 
formulae  for  the  first  order  approximation  by  comparing  the 
terms  of  the  ray  series  with  exact  formulae  for  displacement 
for  reflected  and  refracted  waves  in  a  homogeneous  plane 
parallel  layered  medium.  Unfortunately,  for  inhomogeneous 
media,  such  expressions  are  not  available  and  one  must 
resort  to  numerical  techniques.  However,  Pod'  '  yapol '  ski y  did 
demonstrate  that  the  zeroth  order  approximation  may  not  be 
adequate  for  two  special  conditions.  One  of  these  is  for 
rays  which  have  only  travelled  a  few  wavelengths  from  the 
source.  The  other  is  for  rays  reflected  at  angles 
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tt / 2  with  the  interface  normal,  commonly  called  grazing 
rays.  Woodhouse  (1974),  in  his  doctoral  dissertation, 
computed  the  magnitude  of  the  transverse  component  of 
compress ional  waves  propagating  in  a  radially  symmetric 
earth.  This  work  provides  a  lucid  mathematical  account  of 
the  ray  series  method  although  extensive  numerical  results 
were  absent.  Woodhouse  paid  attention  to  the  additional 
component  of  P  waves,  which  is  a  rather  striking  feature  of 
asymptotic  ray  theory.  From  the  practical  point  of  view, 
however,  the  principal  component  is  more  significant  as  it 
carries  more  energy.  Its  evaluation  must  be  carried  out 
numerically.  To  simplify  the  task,  the  results  will  be 
computed  for  plane  layered  media  with  linear  velocity 
gradients  within  the  layers.  This  reduction  leads  to  simple 
ray  geometries  (circular  arcs)  and  allows  considerable  use 
of  analytic  expressions. 

2 . 2  Basic  equations  of  asymptot i c  ray  theory 

The  basics  equations  of  asymptotic  ray  theory  have  been 
discussed  in  numerous  works.  For  details,  interested  readers 
are  referred  to  Babich  and  Alexseev  (1958),  Alexseev  and 
Gel'chinskiy  (1959),  Alexseev  et  al  (1961),  Kline  and  Kay 
(1965),  Cerveny  and  Ravindra  (1971),  Hron  and  Kanasewich 
(1971),  Woodhouse  (1974),  Cerveny  et  al  (1977)  and  Cerveny 
and  Hron  (1979)  among  others.  We  shall  briefly  derive  the 
equations  for  the  case  whe^e 


' 

' 


K 
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X  /  ju  /  p  and  all  their  derivatives  are  continuous 
functions  of  coordinates.  The  case  where  interfaces  are 
present  will  be  discussed  later  on. 


2.2.1  A  medium  wi th  conti nuous  elastic  parameters 


Substitution  of  equation  2.2  into  2.1  yields 

co  ,  O  '  — >*  ->  — ->• 

l  (iw)  :(iw)  N(Wk)  -  id)M(Wk)  +  L(Wk)  =  o 
k=0 


where 


N(wfc)  =  -pwk+(X+,)(wk.VT)VT+,(VT)^wk 


M(W  )  =  (X+M) (V(Wk*VT)+VT(V*Wk) ) 


+  jh(2  (Vt  •  V)  Wk+  (V2x)Wk)+VX  (Wk*Vx) 


+7«x ( VxxWk) +2 (Vjh • Vx ) Wk 
L(W  )  =  (X+*t)  V  (V-Wk)+«V2Wk+VX  (V*Wk) 

(2.6) 

+V^x  ( VxW,  )  +  2  ( V/U  •  V )  W. 

Since  2.5  must  be  true  for  all  w  ,  the  coefficient  of  any 
power  of  w  must  be  zero,  implying 


N  (W  )  -  M(Wq)  =  0  ^(WQ)  =  0 

N(Wk)  -  M(Wk_1)  +  L(Wk_2)  =  0 


(2.7) 


This  is  a  recurrent  system  of  partial  differential  equations 


• 

•  . 


which  allows  us  to  determine  t  (which  is  invariant  for  all 
terms  in  the  series)  and  all  wk  .  If  we  introduce 
w_1  =  w _2  =  0  we  may  write 


N(w  )  -  M(wk_1)  +  L(wk_2)  =  0  (k>0) 


(2.8) 


Eikonal  equations  for  compress ional  or  shear  waves  can 
be  obtained  as  necessary  conditions  for  the  non-zero 
solution  of  the  system  for  k=0.  These  are 


(Vt)  •  (Vx)  =  a 


-2 


a 


=  (  X+2m ) /p 


(2.9) 


or 


(Vt) • ( Vt)  =  g~2,  32  -  mV P  (2.10) 

a  is  the  velocity  of  propagation  of  compressiona 1  waves, 
sometimes  termed  P,  longitudinal,  irrotational  or 
dilatational  waves,  while  3  is  the  velocity  of  propagation 
of  shear  waves,  sometimes  known  as  S,  transverse, 
equi vo luminal  or  distort ional  waves. 

In  inhomogeneous  media  the  el astodynamic  equation 
cannot  be  decoupled  into  two  wave  equations,  each  describing 
the  propagation  of  compressional  and  shear  waves  (Grant  and 
West  1965,  p.  43-47).  In  the  high  frequency  limit,  however, 


_  . 


- 


. 


V  • 


. 
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two  independent  waves  exist  as  long  as  the  elastic 
parameters  and  all  their  derivatives  are  continuous.  It 
should  not  be  inferred  that  all  particle  motion  associated 
with  the  compressional  wave  is  purely  longitudinal  or 
parallel  to  the  ray.  The  results  of  this  chapter  will 
demonstrate  that  for  the  term  k= 1  an  albeit  small  transverse 
component  normal  to  the  ray  is  present. 

The  eikonal  equations,  2.9  and  2.10  plus  Snell's  Law  at 
boundaries,  provide  us  with  all  the  information  necessary  to 
determine  the  kinematic  properties  of  the  waves,  such  as 
wave  front  geometry,  ray  paths  and  travel  times.  Various 
schemes  have  been  designed  to  solve  these  equations.  Wesson 
(1970,1971),  Cerveny  et  al  (1974),  Popov  and  Psencik  (1978), 
Sorrells  et  al  (1971),  and  Cerveny  and  Hron  (1979)  are  only 
a  few  of  the  major  contributors. 

For  evaluation  of  the  coefficients  of  the  ray  series, 
2.2,  we  introduce  three  orthogonal  right-handed  unit  vectors 
t.  ,  ( i = 1 , 2 , 3 )  which  correspond  to  the  tangent,  binormal 
and  normal  to  the  ray  (see  Figure  2.1).  All  calculations 
performed  in  this  chapter  will  assume  that  the  ray  is  a 
plane  curve  confined  to  the  t^-t^  plane  and  that  the 
vector  t2  points  out  of  the  plane  of  propagation.  It  is 
also  assumed  that  the  z-coordinate  points  downward  into  the 


earth  and  that 


, 


q  •  i 
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is  oriented  such  that  its  projection  onto  the  z  axis 
is  always  positive.  In  this  fashion  all  ambiguity  regarding 
unit  vector  orientation  is  removed.  We  shall  employ  one 
property  of  an  isotropic  medium,  namely  that  is 
parallel  to  . 


We  can  wr  i  te  w, 

k 


in  component  form  as 


(2.11) 


1  ~ 

For  compress ional  waves,  the  component  ti  is  termed  the 

2  3 

principal  component  while  w  t„  +  t-  is  termed  the 
additional  component. 


From  2.6,  2.9,  and  2.11  we  obtain 


N(Wk)-t2  =  (-P+HO.  2)W2 

N(Wk)*t2  =  (-p+iia  2)W2 

Substituting  from  2.8  for  N(w^)  yields 

wk  =  -«2/<X+>I)'(M(Wk_1)  -  L(Wk_2))-t2 

Wk  =  -a2/(A+^)-(M(Wk_1)  -  L(Wk_2))-t3 


(2.12) 


(2.13) 


With  the  geometry  we  have  chosen, 


■ 


■ 
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2 

w.  =0.  Thus  the  additional  component  of  P  waves  is  given 

.K. 

by  t3  . 

1  ~ 

The  evaluation  of  the  principal  component  ti  is 
somewhat  more  involved.  Taking  the  scalar  product  of 
equation  2.7  with  Vi  yields 


M (aW^Vx ) • Vi  =  (L(Wk_x)  -  M (W^  t3) ) 


Vt 


(2.14) 


The  left  hand  side  of  this  equation  can  be  recast  by  using 
some  vector  identities  as 

M(aW  vt)  *Vt  =  2p  (  Vt  •  V )  aW,  +  aW,  (  p  V  t+  (Vt*V)p) 

k  K  K 


Introducing  the  variable  t  as  time  along  the  ray, 


( V  t • V ) aW^  = 


x  d(aW^) 


a 


dt 


(VT,V)P  =  -T  It 


and  hence 


a 


.  ~  d (aVC)  i .  9  i  dn 

M (aW^Vx)  -Vt  =  +  aW^:  pV  x  +  ^ 


a  ~  "  \  a 

Substituting  this  expression  into  2.14  gives  a  linear 


ordinary  differential  equation  for  «w 


1 


d (aW^) 

dt 


+ 


1  dp 
2p  dt 


2 

a 

2  p 


-  M  (W?  t3)  )  -Vt 


(2.15) 
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This  is  one  form  of  the  higher  order  transport  equation. 
Eauation  2.15  enables  us  to  find  recursively  all  values  of 
given  W_.  ,  j  =  0 ,  1 . K- 1  and  T 


An  explicit  solution  of  2.15  can  be  found  by  standard 
techniques  (Simmons  1973).  Denoting  as  some  arbitrary 
initial  point  along  the  ray, 


\h 


W^(t)  =  W^(t  ) 
kv  '  k  '  '~0; 


«/p(V' 

p(t> 


a(to) 

a  (t) 


SU0)  +  2a  ( t ) 


where 


a  (5) 

p(S) 


(  (L(W 


k-1 


0 


\  x 


PL (£T  2 
P  ( t) 


(2.16) 


S(t’) 


exp 

.  \ 

\ 

\ 


1 

2 


a 


(C) V2t 


t 


i 


Let  us  discuss  the  term  S.  For  this  purpose  we 
introduce  ray  coordinates  (q  ,q2,t)  where  ^  ,q2 
specify  the  ray  under  consideration  and  t  gives  the  travel 
time  to  any  point  on  the  ray.  For  a  point  source,  qi'q2 
can  be  regarded  as  take-off  angles  at  the  source.  For  given 
ray  coordinates,  we  can  write  the  location  of  any  point 
along  the  ray  as  x  =  x(q1,q2't)  .  Then  for  a  fixed  t, 


variation  of 


— 


. 

- 
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ql'q2  gives  the  equation  of  a  wave  front.  A  useful 
concept  to  introduce  now  is  that  of  the  ray  tube,  i .e.  the 
family  of  rays,  the  parameters  of  which  are  within  the 
limits  q^'q^  +  dqi;  q2'q2  +  dq2  •  The  infinitesimal  cross 
sectional  area  of  the  ray  tube  can  be  evaluated  as 

dcj(t)  =  |  J  (t)  |  dq1dq2  "here 


and 


J(t) 


V2  T 


9x  x  dx 
3ql  3q2 


_1_  d  ( J/ct) 
Ja  dt 


so  that 


s  (t 1 ) 


a  ( t) 
a  (t 1  ) 


da  ( t 1  )  \  2 
da(t)  j 


An  expression  for  the  principal  component  of  P  waves  from 
equation  2.16  is 


w£(t) 


x  /  a (tQ)  p (tQ)  da (tQ)  % 

Wk^t0)  a (t)  p (t)  da (t) 


+ 


1 


t 

r 


Yu 


2  (p(t)a(t))2  J  \ 

0 


d a(£)  2  a  (£) 
da ( t)  P  (  C) 


(  (L  (W  i)  -  M(W^  t3)).Vi)  d? 


(2.17) 


On  the  basis  of  equation  2.17  we  can  compute  the  principal 
component  at  any  point  t  of  the  ray  if  we  Know  its 
value  at  an  earlier  time  tQ  .  Thus  according  to  asymptotic 
ray  theory,  the  principal  components  of  displacement  at  any 
point  along  the  ray  are  expressed  as  relative  displacements; 
i.e.  relative  to  some  initial  value.  For  our  purposes,  we 


■ : 
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will  give  all  displacements  relative  to  a  displacement  of 
1.0  as  measured  on  the  unit  sphere  surrounding  the  seismic 
source.  Due  to  the  difficulty  of  studying  sources  in 
inhomogeneous  media,  we  further  assume  the  medium  contained 
within  this  sphere  to  be  homogeneous.  To  be  completely 
rigorous,  a  function  describing  the  directivity  or  radiation 
pattern  of  the  source  should  be  included  in  all  expressions 
for  displacement.  However,  the  directivity  will  be  taken  as 
equal  in  all  directions  outward  from  the  source. 


For  the  leading  term  of  the  ray  series,  W 


0 


the 


integral  in  2.17  vanishes  as  does  w 


0 


This  means  that 


_  i  ~  x  '  “(to’  P(V  dCT(V  ~ 

W0  0  1  ’VV  ot  ( t)  p(t)  da  (t ) 


This  zeroth  order  solution  is  exactly  that  obtained  by 
application  of  geometrical  ray  theory.  If  it  is  assumed  that 
energy  propagates  only  along  rays  (the  validity  of  this 
assumption  is  an  entirely  different  topic),  then,  in  the 
zeroth  order  approximation,  no  energy  flows  through  the  side 
walls  of  the  ray  tube.  This  confinement  of  energy  flux  does 
not  hold  in  general  for  any  k,  but  only  for  k=0. 


Similar  expressions  can  be  derived  for  principal  and 
additional  components  of  S  waves.  The  additional  component 
of  an  S  wave  i s 


( X+u ) 


(M(Wk_x) 


(2.18) 


,  i  - 


' 
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while  the  principal  component  of  SV-polarized  waves  is 


3  (tn  )  p(tQ)  da(tn) 


0 


'0 


Wk(t)  “  Wk(t0)  3 ( t )  P(t)  da ( t) 


+ 


+ 


1 


2  ( p  (t)  3  (t)  )' 


'  ZIP'"'1  daja_'2 

\  p (?)  da (t) 


0 

(L(Wk_x)  -  M(Wk)  )'t1  d5 


(2.19) 


It  should  be  pointed  out  that  this  expression  for  the 
principal  component  of  shear  waves  is  valid  only  for  the 
specific  two-dimensional  geometry  employed  herein.  In 
general,  the  principal  component  has  two  vector  components 
and  the  coupling  between  them  may  be  problematic  (Cerveny 
and  Hron  1979) . 


2.2.2  The  effect  of  an  interface 


The  formulae  derived  in  section  2.2.1  can  only  be  used 
when  A,  m  and  p  and  all  their  derivatives  are  continuous 
along  the  ray.  When  the  ray  strikes  a  surface  where  the  k-th 
derivatives  of  the  elastic  parameters  become  discontinuous, 

w.  also  becomes  discontinuous.  It  is  assumed  that  the 

k 

media  on  either  side  of  the  interface  are  in  welded  contact. 
Reflected  and  refracted  compress ional  and  shear  waves  must 
be  introduced  in  order  to  satisfy  the  condition  of  welded 
contact.  The  properties  of  these  waves  can  be  determined 
from  the  natures  of  the  media,  interface  and  incident  wave. 
We  say  that  an  interface  is  of  order  n+1  when  at  least  one 
elastic  parameter  and  its  derivatives  of  order  less  than  n 


- 


.  n  ' 


' 
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are  continuous  across  the  interface,  but  the  n-th  derivative 
i  s  di  scont  i  nuous  (  >  0  )  . 

The  problem  of  reflection  and  refraction  of  elastic 
waves  at  an  interface  is  a  fundamental  one  in  seismology.  A 
host  of  workers  have  devoted  much  time  to  this  problem. 
Several  good  sources  of  information  include  Knott  (1899), 
Zoeppritz  (1919)  and  Ewing  et  al  (1957). 

For  P-SV  waves  propagating  along  plane  rays,  there  are 
four  boundary  conditions  expressing  continuity  of 
displacement  and  stress  across  the  interface.  It  turns  out 
that  for  the  type  of  medium  we  have  chosen,  that  upon 
solving  the  boundary  conditions,  the  case  for  SH  waves  can 
be  treated  independently.  We  shall  employ  this  result  a 
priori  and  omit  further  mention  of  SH  waves.  The  incident 
wave  plus  the  four  waves  arising  at  the  boundary,  reflected 
and  refracted  compress ional  and  shear  waves,  must  satisfy 
these  conditions.  We  shall  designate  these  waves  by  an  index 
n,  where  n= 1  and  2  refer  to  compress i ona 1  and  shear  waves  in 
the  upper  medium  while  n=3  and  4  refer  to  compressiona 1  and 
shear  waves  in  the  lower  medium.  The  incident  wave  is 
labelled  n=0.  Furthermore,  we  assume  that  displacement  for 
each  of  these  waves  can  be  expressed  in  terms  of  a  ray 
series  expansion, 

00  _i 

Wn  =  exp  ( ioo  (t  -  t  )  )  T  W  (i^) 

n  k=0 


(2.20) 


. 


•t  ,*ttf  r 
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where  the  superscript  value  n  refers  to  wave  type,  not 
component.  vn  shall  be  used  to  denote  the  velocity  of  the 
n-th  wave: 


V 


1 


v2=Bi 


Vs 


2 


The  phase  functions  are  given  by  solutions  to  the 

eiKonal  equations  2.9  and  2.10 


(V 


T  ) • CVT 
n 


n=l , 2 , 3 , 4 


(2.21  ) 


with  initial  conditions  on  the  interface  Tn=To  •  We 
employ  local  Cartesian  coordinates  (x,z)  with  positive  x 
pointing  to  the  right  and  positive  z  pointing  into  the 
earth.  The  point  0,  where  (x,z)=(0,0),  is  the  point  where 
the  incident  ray  intersects  the  interface.  Then 


9t  0Tn  sinQ 

n  _  _ 0  _  _ 

9x  9x  Vn 


(2.22) 


where  0n  is  the  acute  angle  between  the  z  axis  and  the  ray 
of  the  n-th  wave.  Equation  2.22  is  merely  Snell's  Law. 


For  waves  leaving  the  interface  into  the  first  medium, 
9t^/9 z<0  for  n=1,2  at  0.  Simi 1 ar ly , 9x^/9 z>0  for 

n=3,4  at  0.  Then  at  point  0, 

9t  e  cos6 

-jf  =  (-1)  n  -y— 51  (2.23) 


where 


. 


■ 

-  *J> 


. 
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e  =1  for  n=l , 2 
n 


If  3Tn/3x>V~1 
u  n 

decrees  that 
introduce 
incident  wave 
respectively. 


e  =  2  for  n=3 , 4 
n 

,  then  sin0^>l  and  the  radiation  condition 

X  X 

2  J- 

co s 0  =  -i(sin  9  -1)  2  .  We  shall  also 

n  n 

where  e0  =1  or  2  depending  on  whether  the 
impinges  the  interface  from  below  or  above, 
Then , 


3t  e  cos0 

u  =  (_1}  u  u 


9  z 


V 


0 


At  the  point  0,  we  introduce  two  perpendicular  unit 
vectors  and  nigV  for  each  n  connected  with  the  rays 

of  the  waves  under  consideration.  m11  is  the  unit  vector 

P 

parallel  to  the  ray  and  pointing  in  the  direction  of 

propagation.  is  perpendicular  to  and  is  oriented 

so  that  its  z  component  is  always  positive.  and  m^7 

correspond  to  the  vectors  and  £3  respectively.  Thus, 


,7n  ,  T7n 

W  m  +  w  m 

k,P  P  k,SV  SV 


(2.24) 


Eight  components  now  exist  for  the  four  waves  generated  at 
the  boundary.  Four  are  principal  components, 


(w 


l 


w; 


w; 


w. 


'k , P  "k,SV  "k,P  '  k,SV) 

and  four  are  additional  components, 

.1  ,,2 


(w: 


w“  wJ  w4  ) 

k,p  k,sv  k,P 


k ,  SV 

Provided  all  lower  terms  in  the  series  are  Known,  the 
additional  components  can  be  computed  at  any  point  along  the 
ray  without  Knowledge  of  initial  values  (see  2.13  and  2.18). 


■ 
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It  follows  then  that  the  only  undetermined  quantities  are 
the  four  principal  components.  These  values  will  serve  as 
initial  values  in  equations  2.17  and  2.19. 


The  projections  of  wn  and  on  the  x  and  z  axes 

will  be  required  for  n=0,1,2,3,4.  We  ascribe  the  variables 

n  r^ri  n 

X  '  wz  and  'flkx>  wkz  to  these  quantities.  After  some 
geometrical  considerations  we  may  write 


6+1 

WV  =  Dsin0  +  (-1)  n  cos0 

kx  k ,  P  n  k ,  SV 

e 

wf  =  (-1)  n  wf  cose  +  W  sine 
kz  k,P  n  k,SV  n 


n 


(2.25) 


The  assumption  that  the  two  solid  elastic  media  are  in 
welded  contact  at  the  boundary  implies  that  stress 
(\z'  nzz^  and  displacement  components  (w x,  wz)  be 

continuous  across  the  interface.  The  relationships  between 
stresses  and  displacements  are 

3W^  3W  ' 

nxz  =M  Tz~  +  ~3x~ 


3W  X3W 


-  (A+2“>  Tz 


X 


3x 


The  four  interface  conditions  can  be  written  as 

l  (-1)  n  wn  =  (-1)  °w° 

n=l 


l  (-1)  n  Wn  =  (-1)  °w° 

u ,  z  z 

n=l 


4  e 

I  (-D  n2Z(w  ) 

n=  1 


(-D  \2(S°) 


. 
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I  (-1) 

n=l 


n  (Wn) 

xz  ' 


(-1) 


'0 


n  (w° ) 

xz  ' 


(2.26) 


We  substitute  in  the  values  of  the  principal  and  additional 
components  from  2.20  and  2.25.  After  a  voluminuous  amount  of 
tedious  algebra,  we  arrive  at  a  set  of  four  simultaneous 
linear  equations  for  the  four  principal  components, 
expressible  as 


wi  th 

A  = 


A  W  =  F 


(2.27) 


sinG 


COS0 


-sinG 


cosG 


-cosG 


sinG 


-cosG 


-sinG 


-plalcos2G  2  P-|$1sin2G2 


p2a2COs264  p2$2S^n2Q  4 


P-l31cos202  Y2p2G2sin2e3  -p232cos2^4 


and 


w1 


w 


w  = 


w 


k  ,P 
2 

k ,  SV 
3 


W 


k,P 

4 

k ,  SV  _ 


and 
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where 


1 


F  = 


F1  =  <_1)  0  WL  -  Wk,SVCOs6l  -  Wk,PSin02 


'Wk,SVCOS03  +  Wk,PSln04 


e  +1 

F2  =  (-1)  0  W°z  -  wliSvsin91  +  W^pcose2 


+  Wk,SVsin03  +  Wk,PCOS04 


eo  ( 

F,  =  (-1)  i  A 


WV  +  2/1  ^ 

3-eQ  3x  kx  y  3-eQ  3-eQ  /  kz 


.  ®0|\  3wk-l,x  ,  ^’k-l.z 

+  ‘-1*  \X3-e0  - ~  +  VA3-e0+2w3-e0. - ~ 

P  / 

-piYiBiwk,svsin26i  +  17^1  -  20isin2e2  w2jP 


p2T282Wk,SVsln203  _  S2  a2  "  262Sln  ®4  ”k,P 

4  e  +1  '  9wf  ,  \  3Wp  _ 
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Note  that  the  matrix  X  is  independent  of  the  value  of 
K.  In  general  to  solve  2.27  for  any  k=j,  it  must  have  been 
previously  solved  for  all  values  k=0 , 1 , . . . , j- 1 .  Only  then 
can  the  variables  F^  be  computed.  It  should  be  stated  that 
for  k=0,  equation  2.27  degenerates  to  the  usual  set  of  plane 
wave  reflection  and  refraction  coefficients  as  discussed  and 
employed  by  numerous  researchers. 

2 . 3  Numer i ca 1  eva 1 uat ion  of  the  term  k= 1  for  £  waves  i n 
medi a  wi th  1 i near  ve loci  tv  gradients 

In  this  section  we  shall  discuss  techniques  of 
numerical  evaluation  of  the  equations  for  principal  and 
additional  components  of  seismic  waves  (2.13,  2.17-2.19)  as 
well  as  evaluation  of  the  boundary  conditions  (2.27).  This 
analysis  will  only  be  performed  for  compress i ona 1  waves. 
More  involved  relationships  can  be  developed  for  shear  waves 
and  shall  thus  be  omitted.  The  key  points  in  need  of  inquiry 
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are  selection  of  a  medium,  rewriting  of  the  vector  operators 
into  a  form  better  suited  to  numerical  computation  and 
evaluation  of  terms  appearing  in  the  boundary  conditions. 

2.3.1  Selection  of  a  med i urn 

One  model  which  lends  itself  to  easy  calculation  is  a 
vertically  inhomogeneous  medium.  If  we  specify  that  the 
velocity  of  P  waves  be  given  by 

a(z)  =  aQ+Va*z,  aQ=constant,  Va=constant 

where  z  is  a  unit  vector  in  the  positive  z  direction,  then 
ray  paths  are  circular  arcs  (Nettleton  1940).  This 
approximation  produces  a  sizeable  reduction  in  the  numerical 
calculations  since  many  of  the  analytical  and  .geometrical 
properties  of  circular  ray  paths  can  be  utilized.  One  such 
saving  is  the  abridged  determination  of  the  ray  take-off 
angle  at  the  source  which  causes  the  ray  to  arrive  at  a 
prescribed  epicentral  distance.  We  shall  also  impose  the 
condition  that  all  geological  interfaces  be  plane  and 
parallel  to  one  another. 

We  shall  utilize  common  relations  for  the  other  elastic 
parameters  as  discussed  by  Gardner  et  al  (1974),  viz., 

3  (x)  =  a  ( z )  /3  2 

P(z)  =  c^(z)3% 


A (z)  =  Z)  =  a9/4 (z)/3^ 
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The  use  of  compact  analytic  expressions  for  the  elastic 
parameters  allows  very  easy  differentiation  which  is 
essential  in  evaluating  the  vector  operators  L  and  M 

The  choice  of  such  models  also  makes  it  possible  to 
write  simple  analytic  expressions  for  the  geometrical 
spreading  of  the  wave  front  dcr  ( tQ) /da  (t)  ( Cerveny  and 
Ravindra  1971,  pp.  86-88). 

2.3.2  The  vector  operators 


For  k=1  the  vector  operators  in  equations  2.13  and  2.17 
must  be  in  a  computationally  viable  form.  Consider  first  of 
all  2.13,  the  additional  component  of  the  P  wave 

W2  =  — (M ( WQ ) *t3)a2/(A+u) 


The  component  of 
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For  P  waves  Vt  is  parallel  to  t  and  we  can  state 

V  t  =  t^/'a 

Expansion  of  equation  2.17  and  using  2.6  and  2.28  shows  that 
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where  K  is  the  ray  path  curvature.  K  is  taken  as  positive  if 
the  center  of  the  circle  describing  the  ray  at  a  point  lies 
stratigraphical ly  above  the  ray. 


The  derivatives  of  Q  where  Q  is  any  of  the  elastic 
parameters  a,p,A,y  can  be  performed  easily  with  the 
medium  discussed  in  the  previous  section.  Let  T  be  the  angle 
between  the  positive  z  axis  and  the  unit  vector 

£]_(0-T-tt)  .  Then  it  holds  that 


dQ 
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and  the  derivatives  with  respect  to  z  are  readily  available. 


The  t1  and  derivatives  of  wj  require  more 
care.  Once  the  ray  has  impinged  upon  a  boundary,  it  becomes 
impossible  to  write  an  expression  for  the  displacement  which 
lends  itself  to  analytic  differentiation.  For  this  purpose, 
we  propose  the  following  algorithm: 

(1)  Determine  the  ray  take-off  angle  at  the  source 
which  causes  the  ray  to  emerge  at  a  prescribed 
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epicentral  distance.  We  term  this  ray  the 
central  ray. 

(2)  Partition  each  central  ray  segment  into  NPT+1 
equal  travel -time  subsegments  thus  defining 
NPT  points  along  each  segment. 

(3)  Shoot  two  rays  each  to  the  left  and  right  of 

the  central  ray.  Locate  the  point  where  each 
of  these  auxiliary  rays  intersect  the  plane 
normal  to  the  central  ray  at  all  NPT  points 
(see  Appendix  A  for  algebraic  details).  At 
these  points  the  displacement  can  be 

computed  and  by  fitting  a  spline  through  the 
five  ( 1 3 ^ wq ( 1 3 ) )  points,  values  for 

3WQ/3t  and  can  be  calculated  (see  Figure 

2.2)  . 

1  2  1  2 

(4)  For  computation  of  and  3  wr/3tl 

at  the  j-th  point  of  the  ray  segment,  employ 
values  of  wj  at  the  j-2,  j-1,  j,  j+1,  j+2-th 
points  of  the  ray  segment.  Consult  Figure  2.3 
for  details.  A  spline  can  be  fitted  through 
the  five  points 

t,  ,  wj 

1  ,m  0  ,m 

and  the  derivatives  at  m=j  determined.  Of 
course,  this  procedure  requires  an  obvious 
modification  if  j<2  or  j>NPT-2. 

Note  that  the  methods  outlined  for  computing  directional 
derivatives  of  displacement  will  only  give  good  results  for 
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Figure  2.2  Evaluation  of  the  t.,  derivatives. 


34 


small  values  of  ray  curvature.  For  this  reason,  we  have 
limited  our  study  to  rather  weak  inhomogeneities.  If  one 
were  to  extend  these  approaches  to  strongly  inhomogeneous 
media,  the  technique  would  require  the  ray  segmentation 
(step  2  above)  to  be  finer  for  greater  degrees  of 
i nhomogenei ty . 

Let  tQ  be  the  initial  time  point  of  the  ray  segment 
or  the  time  to  reach  the  unit  sphere  surrounding  the  source 
if  the  very  first  ray  segment  is  under  consideration,  and  t 
be  the  time  at  which  the  central  ray  segment  terminates  on 

an  interface.  Then  to  evaluate  the  integral  in  equation  2.17 

->  ->  ->  ->■  , 

from  tQ  to  t ,  the  vector  operators  L-t^/a  and  M»t^/a 
should  be  computed  at  the  NPT  points,  and  a  spline  passed 
through  them.  In  this  manner,  the  vector  operators  can  be 
determined  by  substitution  for  any  time  £,tQ^<t  ,  along 

the  ray  segment.  Since  the  quantities  da(£)/da(t), 
a(£)r  p(£)  can  be  given  by  analytic  expressions, 

integration  of  2.19  is  simplified.  By  computing  values  of 

1  3 

wk  and  on  the  incident  side  of  interface  we  have  the 

necessary  requirement  to  solve  equation  2.27. 

2.3.3  Evaluation  of  boundary  condi tions 

Matrix  A  in  equation  2.27  is  dependent  only  upon  the 
elastic  parameters  on  both  sides  of  the  interface  and  upon 
the  angle  of  incidence.  In  the  variables 
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F1  to  F4  it  is  necessary  to  Know  the  values  of  the  four 
additional  components  corresponding  to  reflected  and 
refracted  waves  at  the  interface.  These  can  be  computed  by 
equations  2.13  for  P  waves  (superscript  n= 1,3)  and  2.14  for 
S  waves  (superscript  n=2,4). 

The  evaluation  of  the  x-der i vat i ves  of 
( n=0 , 1 , 2 , 3 , 4 )  are  also  facile.  Since  we  actually  have  five 
rays  impinging  on  the  interface,  the  boundary  conditions  for 
k=0  must  be  solved  five  times.  The  x-derivative  problem  then 
is  a  simple  one  of  spline  fitting  through  the  five 
solutions . 

The  z-derivatives  of  ( n=0 , 1 , 2 , 3 , 4 )  are  not  quite 

as  accessible.  According  to  the  zeroth  order  approximation, 
as  a  wave  moves  away  from  the  interface,  the  only  factor 
affecting  w  is 

/CL  .  P  .  ,  ,  .  \  \  .  • 

/  1  1  da (u) 

a(z)  p(z)  da (z) 

where  the  subscript  " i "  denotes  the  value  of  the  variable  at 
the  interface  and  "u"  represents  the  unit  sphere  surrounding 
the  source.  " RC"  is  the  angle  dependent  reflection 
coefficient,  available  from  the  solution  of  equation  2.27 
for  k=0.  For  the  wave  travelling  away  from  the  interface,  it 
is  necessary  to  compute  an  expression  relating  the 
displacement  of  the  wave  at  a  vertical  distance  z  from  the 
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interface  to  the  displacement  at  the  interface,  z=0. 
Following  Cerveny  and  Ravindra  (1971,  pp.  74-88)  for  our 
medium  and  omitting  the  subscript  0  and  superscript  n 


W  ( z )  =  W 


(v^p^sinG  ) 

1  -s  ^4  /  / 

3  (v  (z)  ) 


r  ( z)  drA~--  (l-v2(z)SV2)  1 


99 


0 


(2.31 ) 


where 

w. 

i 


e 


o 


v 


0 

K 


SV 


.  displacement  at  the  boundary  after  reflection  or 
ref ract ion , 

.  velocity  and  density  experienced  by  the 
incident,  reflected,  or  refracted  wave  at  the 
boundary , 

.  ray  take-off  angle  with  the  positive  z  axis 
measured  at  the  source 

.  velocity  at  the  source, 

.  velocity  gradient  in  the  layer  in  which  the  wave 
propagates , 

z 

sin0o/vo  r(z)  =  +  j  tanG(z)dZ 

0 


v  ( z)  = 


v  .  +K: 

i 


tane(z)  =  v(z) (1-v2 (z) SV2)  "  SV 


...  horizontal  distance  travelled  from  source  to  point 


of  incidence  with  the  interface. 

9W  (z 

The  problem  then  is  to  evaluate  D  =  --dz— 


z=0 
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To  accomplish  this,  one  can  find  the  x  distance  where  the 
shallowest  propagating  ray  impinges  the  interface,  say  at 

X  =X0  2.32 

The  z  coordinates  where  the  other  (four)  rays  intersect  this 
vertical  line  are  determined  and  the  zeroth  order 
displacements  calculated  there.  The  desired  z  derivative  is 
then  obtained  by  spline  fitting. 

It  should  be  indicated  that  as  the  angle  of  incidence 
of  waves  upon  the  interface  approaches  a  critical  angle,  the 
values  of  the  z-der i vat i ves  of  critically  refracted  waves 
tend  to  infinity.  For  this  reason  only  subcri tical ly 
reflected  rays  will  be  considered  in  our  numerical 
calculations.  This  singularity  will  clearly  be  observable  in 
the  results  of  the  next  section. 


2.4  Numerical  results 


In  a  medium  composed  of  a  laterally  homogeneous  layer 
overlying  a  laterally  homogeneous  halfspace,  the  equations 
discussed  in  the  previous  section  can  be  evaluated 
numerically  to  give  an  indication  of  the  magnitudes  of  the 
terms  k=0 , 1  in  the  ray  series.  For  each  of  several  models 
studied,  four  graphs  are  shown.  The  first  is  a  plot  of  the 
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magnitude  of  wQ  =  versus  distance  from  the  source.  Also 

1  3 

given  are  plots  of  the  magnitudes  of  and  w  versus 
distance  from  the  source.  In  all  cases,  the  displacements 
are  relative  to  a  unit  displacement  on  the  unit  sphere 
surrounding  the  source.  The  fourth  graph  is  a  plot  of  U= 

|  w 1 / ( 2tt  ! wQ  |  )  x  ioo%  .  This  is  a  ratio  between  the 
magnitudes  of  the  K= 1  and  k=0  terms  in  the  ray  series.  The 
greater  this  value,  the  greater  the  effect  of  the 
non-leading  term  in  the  ray  series.  This  ratio  must  be 
divided  by  the  source  frequency,  f,  in  Hertz,  to  obtain  the 
frequency  dependence.  Thus,  reference  to  a  value  of  U  such 
as  25/ f% ,  for  example,  means  that  for  source  frequency  of, 
say, 10  Hertz,  the  ratio  of  the  k=1  term  to  the  k=0  term  is 
25/10=2.5%.  Also  given  in  each  set  of  graphs  are  velocity 
depth  profiles  of  the  media  under  consideration.  The  label 
numbers  on  these  profiles  correspond  to  those  on  the 
individual  curves.  Note  that  all  reflections  are  subcritical 
in  keeping  with  the  discussion  of  the  previous  section. 

In  the  first  model,  we  study  P  waves  propagating 
through  a  homogeneous  layer  and  reflecting  from  a 
homogeneous  ha  If space.  The  refractive  index  at  the 
reflecting  interface  is  assigned  an  increasing  value  in  each 
of  the  four  cases  with  the  interface  becoming  nearly 
transparent  in  the  last  one.  These  results  are  displayed  in 
Figure  2.4.  The  curves  for 
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j w q  |  behave  as  expected  with  the  amount  of  reflected  energy 
being  directly  linked  with  the  refractive  index  across  the 
reflector.  As  well,  the  derivatives  of  the  |wj|  curves  tend 
to  infinite  values  as  one  approaches  the  critical  distance. 
It  is  interesting  to  examine  the  curves  for  |w^j  in 

conjunction  with  those  for  the  additional  component  |w^| 
noting  that  the  latter  varies  roughly  as  the  derivative  of 
the  former.  In  general,  the  points  where  the  curves  for 

1 1  |  3 1  i 

I  and  |W1|  dip  to  0  indicate  approximate  points  where  wl 

3 

and  have  polarity  reversals.  The  curves  for  the  ratio  U 
indicate  that  for  rays  incident  on  the  reflector  at 
substantially  less  than  critical  angles,  the  magnitude  of 
the  k=1  term  is  approximately  25/f%  for  all  but  the  most 
transparent  interface. 

The  second  test,  displayed  in  Figure  2.5,  examined  P 
waves  propagating  through  an  inhomogeneous  layer  and 
reflected  from  a  homogeneous  ha  If space.  Both  positive  and 
negative  velocity  gradients  in  the  layer  were  considered.  In 
order  to  retain  some  degree  of  similarity  between  the 
models,  the  refractive  index  at  the  reflecting  interface  was 
kept  constant.  Once  again,  note  the  similarity  between  the 

i  3 i  ill 

curves  for  | |  and  the  derivatives  of  the|WQj  curves. 
However,  it  must  be  recalled  that  in  inhomogeneous  media, 
depends  on  spatial  derivatives  of  the  elastic 

parameters  as  well  as  the  derivatives  of  the  zeroth  order 
approximation  wj  , 
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The  zeroes  in  the  lW]J  ,  I  I  and  ^  correspond  to 
polarity  reversals  in  the  functions.  The  singularity  in  the 
term  K= 1  as  the  receiver  approaches  the  critical  distance  is 
evident.  However,  for  subcritical  reflection  in  the  model, 
the  effect  of  the  non- leading  term  in  the  first  order 
approximation  is  no  more  than  25/f%.  It  can  also  be  inferred 
from  inspection  of  these  curves  that  as  one  increases  the 
value  of  the  absolute  velocity  gradient  (and  hence  the 
degree  of  inhomogeneity  of  the  medium),  the  magnitude  of  the 
k= 1  term  also  increases. 

A  commonly  used  technique  in  refraction  seismology 
employs  refracted  waves  that  occur  when  the  velocity  varies 
continuously  with  depth.  It  is  interesting  to  study  the 
relationship  between  the  geometrical  ray  theory  result,  |wj 
W.,  |  and  |wj|  for  refracted  rays  in  an  inhomogeneous  ha  If  space. 
Computed  results  for  such  an  inquiry  are  shown  in  Figure 

i  1 i  | w3 |  ili 

2.6.  Curves  for  I wq  1  )|VVi'and  | W  J  all  have  the  same  basic 
nature,  namely,  monotonic  decrease  with  distance  from  the 
source.  The  graph  of  the  ratio  U  versus  distance  indicates 
two  conclusions.  First,  that  the  effect  of  the  K= 1  term  is 
rather  small  for  rays  refracted  through  an  inhomogeneous 
half space  and  second,  that  the  ratio  of  magnitudes,  U, 
varies  almost  linearly  with  the  degree  of  inhomogeneity. 

The  zeroth  order  approximation  to  the  solution  of  the 
el astodynamic  equation  predicts  that  the  amplitude  of  a  wave 
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reflected  from  an  inhomogeneous  halfspace  will  not  depend  on 
the  degree  of  inhomogeneity  of  the  ha  If space.  In  the  first 
order  approximation  this  is  no  longer  so.  Brekhovskikh 
(1960)  studied  a  similar  problem  for  plane  electromagnetic 
waves  propagating  through  a  homogeneous  layer  and  reflected 
from  an  inhomogeneous  ha  If space.  For  spherical  elastic  waves 
propagating  in  an  elastic  solid,  no  previous  results  have 
been  published  on  a  similar  problem,  to  the  author's 
know  1  edge . 

The  next  model  tested  considers  P  waves  propagating  in 

a  homogeneous  layer  and  reflected  from  an  inhomogeneous 

ha  If space.  The  degree  of  vertical  inhomogeneity  of  the 

ha  If space  has  been  varied,  but  the  refractive  index  across 

the  interface  has  been  kept  constant.  With  this  condition, 

it  is  not  surprising  to  notice  in  Figure  2.7  that  the  |wj| 

i  3  i 

curve  for  all  seven  cases  is  identical.  |W,  j  ,  since  it 
depends  on  the  derivative  of  wj  ,  behaves  similarly. 

1  w  |  and  the  ratio  between  the  terms,  U,  are  however, 
inhomogeneity  dependent.  With  an  increase  in  the  value  of 
the  velocity  gradient  of  the  halfspace,  both  |w||  and  U 
increase.  The  singular  nature  of  the  curves  as  the  receiver 
approaches  the  critical  point  is  also  evident. 

In  this  model  it  is  instructive  to  examine  the 
solutions  to  the  boundary  conditions  (equation  2.27)  for  k=0 
and  k=1.  These  are  displayed  in  Figure  2.8  for  the  models 
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under  consideration.  Using  the  notation  of  equation  2.27, 
the  curves  give  the  ratio  of  the  reflected  compressiona 1 
displacement,  p  to  the  incident  compress ional 

displacement,  for  k=0,1.  This  ratio  for  k=0  is  more 
commonly  known  as  the  reflection  coefficient.  For  k=1,  this 
simplicity  is  lost  due  to  the  increased  complexity  of  the 
higher  order  boundary  conditions. 

It  is  desirable  to  remove  the  effect  of  the  singularity 
at  the  critical  point  and  be  able  to  study  the  effect  of  the 
k= 1  term  for  rays  which  are  incident  upon  plane  interfaces 
at  angles  near  90°  .  To  achieve  this  the  programs  have  been 
run  for  a  homogeneous  layer  overlying  an  inhomogeneous 
halfspace  with  the  velocity  in  the  layer  exceeding  that  at 
the  top  of  the  halfspace.  These  results  are  presented  in 

Figure  2.9.  Again,  for  all  seven  velocity  depth  profiles, 

i  1  ■  i  3 

the  curves  for  wn  I  w,  I  are  the  same.  The  results  for 
drop  sharply  with  distance  from  the  source,  pass 
through  zero  between  30  and  50  km  and  decay  very  slowly  as 
the  rays  approach  grazing  incidence.  The  effect  of  this  upon 
the  curve  of  the  ratio,  U,  between  the  magnitudes  of  the  k=0 
and  k= 1  terms  is  to  gradually  increase  as  the  grazing  angle 
is  approached.  It  is  evident  that  in  this  case,  the  effect 

of  the  term  w  is  less  than  50/f%. 

1 

In  homogeneous  media,  a  tractable  expression  can  be 
derived  for  the  additional  component  of  once  reflected  P 


-  a 


km/sec 


48 


£_  >>  <D 
O  —  JO 


a> 


<D  <D  aJ 


cr> 


49 


waves.  If  R(0)  is  the  plane  wave  reflection  coefficient  as 
a  function  of  incident  angle,  0  ,  (see  the  solution  of 
equation  2.27  for  k=0)  then  some  elementary  algebra  and 
geometry  shows  that 

2 

3  _  acos  0  dR ( 0 ) 

1  ‘  4h2  d9 

where  a  is  the  P  wave  velocity  in  a  layer  of  thickness  h. 
The  evaluation  of  this  expression  is  shown  by  the  dotted 
line  in  Figure  2.9.  Clearly,  the  matching  between  analytic 
and  numerical  techniques  is  very  good  in  this  simple  medium. 

The  solutions  to  the  boundary  conditions  for  this  model 
are  presented  in  Figure  2.10.  It  can  be  shown  that  the 
zeroth  order  solution  to  the  boundary  conditions  tends  to  -1 


as  the 

incident  angle 

approaches 

CD 

O 

0 

.  The 

first  order 

solution 

increases  without 

limit 

for 

rays 

approaching 

grazing 

angles.  However, 

the 

|Wd| 

curves 

remain 

finite  due 

to  the  fact  that  the  geometrical  spreading  term  da(tQ)/dj(t) 
in  equation  2.17  tends  to  zero  at  a  faster  rate. 


2.5  Conclusions 


A  commonly  used  technique  for  describing  the 
propagation  of  seismic  body  waves  in  an  inhomogeneous  medium 
is  the  asymptotic  ray  series.  Most  seismologists  have  been 
satisfied  with  the  zeroth  order  approximation  which  is 
correct  for  co  -+■  , 
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Figure  2.10  Solutions  to  the  k=0  and  k= 1  boundary  conditions 
for  P  waves  reflected  from  a  vertically  inhomogeneous 
halfspace.  The  velocity  distribution  at  the  interface  is 
given  in  Figure  2.9. 
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This  relies  on  the  assumption  that  the  higher  terms 
in  the  ray  series  (k£1)  could  be  neglected.  This  chapter  has 
presented  the  formulae  for  the  first  order  approximation  in 
a  numerically  practical  form.  These  equations  have  been 
evaluated  for  some  common  problems  in  re  lection  and 
refraction  seismology.  On  the  basis  of  these  tests  we  wish 
to  establish  a  "rule  of  thumb"  on  the  magnitude  of  the  k=1 
term  relative  to  the  leading  term  which  embodies  geometrical 
ray  theory.  For  subcritical  reflection  from  an  interface 
separating  vertically  inhomogeneous  media,  the  error  in 
neglecting  the  k=1  term  is  generally  less  that  75/ f%  where  f 
is  the  source  frequency  in  Hertz.  For  refracted  arrivals  the 
error  is  considerably  less  than  that.  Since  a  reasonable 
seismic  source  frequency  for  crustal  explosion  seismology  is 
of  the  order  of  10  Hertz,  the  error  is  less  than  8%.  For 
exploration  seismic  work,  the  error  drops  even  lower  to 
1-5%,  which  will  be  imperceptible  on  individual  seismic 
records . 

This  problem  of  higher  terms  in  the  ray  series  could  be 
further  researched  by  extending  the  presented  techniques  to 
more  complicated  media.  One  of  these  would  be  to  generalize 
the  processes  to  an  arbitrary  number  of  layers.  Then  the 
first  order  approximation  could  be  studied  for  waves 
propagating  in  layers  that  are  thin  compared  to  the 
wavelength.  The  effect  of  interface  curvature  on  the  first 
order  approximation  could  also  be  investigated.  Results  for 
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the  zeroth  order  approximation  in  media  with  curved 
boundaries  are  presented  in  Chapter  6  of  this  thesis. 
Finally,  the  effect  of  the  k= 1  term  may  play  a  considerable 
role  in  diffraction  phenomena,  since  these  are  known  to  be 
generally  frequency  dependent. 


3.  NUMERICAL  CONFIRMATION  OF  METHODS  OF  COMPUTING  SEISMIC 


WAVE  AMPLITUDE 


3 . 1  Introduction 


A  number  of  geophysical  programs  today  are  employing 
asymptotic  ray  theory  as  a  primary  means  of  theoretically 
investigating  seismic  wave  propagation.  Notable  among  these 
are  Cerveny  and  Ravindra  (1971),  Hron  and  Kanasewich  (1971) 
and  Cerveny  et  al  (1977).  The  reasons  for  the  widespread  use 
are  the  theory's  versatility  in  complex  media  and  the  great 
computational  efficiency  over  full  wave  techniques  such  as 
the  method  of  finite  differences.  It  is  important  to 
determine  the  degree  of  accuracy  of  the  ray  theory  in  order 
to  justify  its  extensive  uses  in  both  crustal  studies  and 
oil  exploration  geophysics.  This  accuracy  may  be  determined 
by  direct  numerical  integration  of  the  displacement 
integral s . 

A  recent  study  (Marks  1976,  Marks  and  Hron  1979a)  has 
presented  a  technique,  based  on  the  parabolic  cylinder 
function  (or  Weber  function)  ,  which  eliminates  the 
singularity  that  asymptotic  ray  theory  produces  when  the 
receiver  lies  at  the  critical  point  of  a  head  wave.  Direct 
numerical  integration  can  also  be  utilized  to  verify  the 
amplitudes  produced  by  the  Weber  function  approach.  At  this 
point  the  reader  may  find  it  beneficial  to  consult  Appendix 
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B  for  a  review  of  the  principles  of  the  Weber  function 
approach  in  pi ane- 1 ayered  media.  Asymptotic  ray  theory  was 
discussed  in  Chapter  2.  We  shall  employ  only  the  zeroth 
order  approximation  here. 


3 . 2  Theory  for  numer i ca 1  i nteqrat i ons 


According  to  Brekhovskikh  (1960)  or  Ewing  et  al  (1957), 
among  others,  the  following  integral  expression  describes 
the  potential  at  the  earth's  free  surface  of  a  P  wave  with 
bottoming  reflection  on  the  L-th  interface  of  a  plane, 
para  1  lei  - 1 ayered  medium  with  all  layers  homogeneous  and 
i sotropic 


-icot 


$ (r , z)  = 


R (aLq/a1) a (aLq/a1) Hq (krq) 


—  2 

exp(ik  B ( q) ) q ( 1-q  )  ‘dq 


(3.1) 


where 

k=to/a 


1‘ 


to  =  2Trf  = 


a  .  /  b  . 
l  i 


R.  = 


a  - 


wave  number  in  first  layer, 
circular  frequency  of  source, 

compress i ona 1  and  shear  wave  group  velocities  in 
the  i-th  layer, 

plane  wave  reflection  coefficient  appropriate  to 
the  type  of  reflection  occuring  at  the  L-th 
interface , 

product  of  all  other  reflection,  transmission,  free 
surface  and  surface  conversion  coefficients 


surface  conversion 


3J  Y1C<y.  3  vsi  J  >Y*A  f  <3  ft 


' 


55 


H0 

B(q) 


encountered  by  the  ray  in  its  propagation  from 
source  to  receiver, 

-  HanKel  function  of  the  first  Kind  and  zeroth  order, 

=  N^h.,  (a2/b2-q2)  2  +  ( (ND1~1) hn +z )  • 


1 


SI  1  1  1 

9  p  r  ,  ,  2  ,  2  2 ,  '2 

(1-q2)  2  +  l  Nn^hi(al/ai“q  } 

L  i=2 

+  J2  Nsihi  (ai/ bi-q2)lj 


=  phase  function  of  the  ray  with  P  type  ray 

segments  and  n  .  S  type  ray  segments  in  the  i-th 

31 

1 ayer , 

=  thickness  of  the  i-th  layer, 
z  =  height  of  receiver  above  the  first  (non-free) 
interface . 


By  application  of  the  well  Known  equations,  one  may  write 
the  expressions  for  the  horizontal  (u)  and  vertical  (w) 
components  of  displacement  as 


u 


-1 


ke 


-icot 
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exp (ikB (q) ) q  l  q 


/ 

1 

1 


where 


56 
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H]_  is  the  HanKel  function  of  first  Kind  and  order  and  the 
arguments  have  been  dropped  from  R  and  a  . 

By  the  proper  choice  of  integration  contours  and  by 
limiting  ourselves  to  high  frequencies,  the  HanKel  functions 
may  be  expanded  asymptotically  as  (Abramowitz  and  Stegun 
1965) 


With  these  approximat ions  and  by  placing  the  receiver  at  the 
free  surface  (  z  =  h-j_  ),  the  displacements  may  be  rewritten 
as 


- —  i  exp  ( -iojt+iTr/4 )  Rcr 
2irr  } 


n 


(3.2) 


where 


B  =  rq 


(3.3) 


The  path 
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D  utilized  in  equation  3.2  must  be  decomposed  depending 
on  the  receiver's  location  relative  to  the  critical  point. 
For  epicentral  distances  less  than  the  critical  point  the 
path  consists  of  two  parts:  path  dq  and  path  D1  (see 
Figure  3.1).  The  former  is  a  semicircle  in  the  upper  half 
plane  of  the  complex  variable  q  ,  while  the  latter  is  given 
parametrically  by 

( l~q2 ) 2  =  (1-x2)  +  p • exp ( -iff/ 4 )  /o  4\ 


for  real  p  in  the  domain  (-00,03)  .  The  variable  x,  which 
is  the  real -valued  saddle  point  for  the  ray,  may  be  obtained 
by  solving  the  condition  of  steepest  descent 


dB 

dq 


=  r 


q=x 


Upon  performing  this  operation,  we  arrive  at  an  equation  for 
rays . 


Beyond  the  critical  distance,  r*  ,  one  must  consider 

2  2  o  J- 

the  branch  cuts  of  the  radical  (a  /v  -  q  )2  which 
correspond  to  the  generation  of  P-type  (if  v=a)  or  S  type 
(if  v=b)  head  waves  from  the  L-th  interface.  In  order  that 
the  deformed  path  D  end  on  the  same  sheet  of  the  Riemann 
surface  as  it  began  on,  it  is  necessary  to  circumvent  the 
cuts  due  to  the  above  radical  by  path 


^  — 


58 


59 


D 


(see  Figure  3.1b).  This  cut  is  chosen  as 


2 


(3.5) 


for  real  values  of  p  in  the  domain  [0,  °°)  and  for 
x*  =  a1/vL+1  ,  Physically,  it  is  clear  that  x*  is  the 
critical  ray  parameter  in  the  first  layer.  The  radical  will 
have  a  different  sign  on  either  bank  of  the  branch  cut  since 
it  crossed  onto  a  different  Riemann  sheet. 

The  use  of  the  contours  of  integration  shown  in  Figure 
3.1  and  described  parametrical ly  by  equations  3.4  and  3.5 
have  been  discussed  at  length  by  many  authors.  Rigourous 
justifications  of  these  paths  can  be  found  in  Cerveny  (1962, 
1965,  1967)  and  Marks  and  Hron  (1979a).  Recently,  Daley  and 
Hron  (1979)  studied  a  similar  problem  for  SH  waves  in 
transversely  is tropic  media. 

For  epicentral  distances  less  than  the  critical  point, 
only  the  reflected  wave  will  exist.  In  this  case,  the 
exponential  term  in  the  integrand  may  be  expanded  as 


+  .  .  . 


where 


,2,2  _ w  2  ,  2  2. -3/2 
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We  also  note  that  b  x)  elated  to  the  two  way  travel 
time  of  the  reflected  ray  by  B(x)/a1  =  tr  .  This  yields  a 

result  which  lends  itself  to  numerical  integration 

1  u  -ia)(t-t  ) 

=  e  r 

w  i 


2Trrx 


R  a 


exp (-k (r-r) p2/2x3) 


\  2  ^ 

(1-q  )2; 


(3.6) 


Along  the  circumvention  of  the  branch  cut  a  similar 

expansion  may  be  produced  as 

ikB  =  ikB  (x*)  -  k  ( 1-x*2)  (r-r*)  ( 1+i)  p/2'2x* 

-k (r-r*)  p2/2x*3  +  ... 

where  r*  is  just  f  evaluated  at  x=x*  .  Again,  b(x*) 

is  related  to  the  travel  time  of  the  head  wave  by 

B  (x* )  /a,  =  t.  • 

1  h 

The  plane  wave  reflection  coefficient  R  can  be  cast  in 
the  form 

R  =  E  -  Dfa^^  -  q2)^  =  E  -  DS  ,g-7) 

where  E  and  D  are  functions  of  a  q/a  .  E  and  D  do  not 

Li  L 

depend  on  the  radical  S,  but  only  on  even  powers  of  S.  By 
substituting  equations  3.5  and  3.7  into  3.2,  and  recalling 
that  the  sign  of  S  differs  from  left  to  right  bank  of  the 
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path,  we  arrive  at  an  integral  expression  for  the  head  wave 
displacement : 


U' 


l  W' 


-iU(t-th)  /  k  \h  r 
-e  - - -  z 


,  2Trrx* 

\ 

\ 


DSa 


0 


.2,1 


exp (-k ( 1-x*  )  2  (r-r*) (l+i)p/2  x* 


r 

q 


k(r-r*)  p  / 2x*  )  9  ^  dP 

i  d-q  )  2 


(3.8) 


Equations  3.6  and  3.8  are  the  final  expressions  that  we 
wish  to  evaluate  along  the  integration  paths  shown  in  Figure 
3.1.  Care  must  be  taken  to  ensure  that  the  correct  sign  is 
used  for  all  radicals  in  the  reflection  and  transmission 
coefficients  in  the  integrands.  The  integrands  are,  in 
general,  complex  valued  and  thus  the  integrals  must  be 
evaluated  once  for  the  real  part  and  once  for  the  imaginary 
par  t . 


3 . 3  Numerical  resu 1 ts 


An  efficient  and  versatile  algorithm  for  numerical 
integration  employs  cautious  Romberg  extrapolation,  ably 
discussed  by  de  Boor  (1971).  On  the  basis  of  a  FORTRAN 
program  with  the  algorithm  at  its  hub,  amplitude  -  distance 
curves  were  obtained  for  various  rays  and  source  frequencies 
for  the  medium  listed  in  Table  3.1.  This  medium  is  basically 
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that  due  to  Camming  and  Kanasewich  (1966).  The  top  two 
layers  have  been  thickened  to  increase  the  ratio  of 
thickness  to  wavelength  to  values  acceptable  by  asymptotic 
ray  theory. 

The  cost  of  computing  results  via  numerical  integration 
exceeds  those  derived  by  asymptotic  ray  theory  or  the  Weber 
function  by  a  factor  of  160  to  200,  depending  on  the 
frequency  and  ray  path  selected.  However,  this  integration 
technique  may  be  taken  as  exact  and  the  Weber  function 
approach  and  asymptotic  ray  theory  will  be  shown  to  be  very 
economical  approximate  methods. 

Figures  3.2  -  3.6  are  a  series  of  plots  of  vertical 
amplitudes  versus  distance  from  the  source  calculates  by 
means  of  asymptotic  ray  theory,  the  Weber  function  approach 
for  interference  reflected  head  waves  or  direct  numerical 
integrations.  Each  diagram  has  a  schematic  sketch  of  the 
ray(s)  under  consideration.  A  solid  line  in  this  sketch 
implies  propagation  in  the  compress ional  (P)  mode,  while  a 
broken  line  denotes  the  shear  (S)  mode.  The  predominant 
source  frequency  (F)  is  indicated  where  necessary. 

Figure  3.2  displays  an  amplitude  distance  curve  for  an 
unconverted  P  wave  reflected  from  the  third  interface.  It 
can  be  observed  that  the  amplitude  curves  derived  by  direct 
numerical  integration  (numbers  2  -  4)  preserve  the  Carrie 
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TABLE  3.1 
The  Alberta  Model 


P  Wave 
Velocity 
km/ sec 

S  Wave 

Velocity 

km/sec 

Volume 

Density 

g/cm3 

Layer 

Thickness 

km 

2.31 

1.  33 

2.04 

8 

3.06 

1.77 

2.21 

8 

6.50 

3.75 

2.73 
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general  shape  of  the  ray  theoretical  curve  (number  1).  At 
higher  frequencies,  the  integration  curves  attain  sharper 
maxima  and  more  closely  follow  the  ray  theoretical  curve. 
The  good  fit  between  phase  shifts  between  curves  1  and  4 
should  also  be  noted. 

In  Figure  3.3,  the  common  logarithm  of  amplitude  is 
displayed  against  the  distance  for  a  reflected  and  converted 
wave.  The  location  of  the  critical  points  for  the  SPS  and 
SSS  head  waves  (marked  Rep  and  Res,  respectively)  are 
indicated.  Following  this,  Figure  3.4  displays  the  amplitude 
of  an  SSS  head  wave,  which  has  its  critical  point  at  the 
location  marked  by  Res  in  Figure  3.3.  The  end  of  the 
interference  zone  between  the  reflected  and  head  waves  is 
indicated  by  Ri .  Thus,  outside  this  region,  there  is  a  near 
perfect  fit  between  the  results  of  the  two  methods. 

Figure  3.5  shows  confirmation  of  the  frequency 
dependence  of  head  waves.  Curves  1  and  2  were  computed  by 
ray  theory  which  predicts  an  inverse  frequency  dependence 
for  head  wave  amplitude.  Curves  3  and  4,  calculated  by 
direct  integration,  produce  a  very  good  fit  with  asymptotic 
ray  theory  outside  the  interference  zones.  The  end  of  the 
interference  zone  for  the  lower  (higher)  frequency  is 
indicated  by  the  rightmost  (leftmost)  Ri .  The  phase  shifts 
are  also  shown  corresponding  to  curves  2  and  4. 
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Figure  3.5  Frequency  dependence  of  head  waves  computed  by 
different  means. 
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Within  the  interference  zone,  it  is  necessary  to 
consider  the  composite  of  the  reflected  wave  and  head  wave 
since  they  do  not  exist  independently  there.  It  was  for  the 
purpose  of  investigating  the  interference  reflected  head 
wave  that  the  Weber  function  approach  was  developed.  It 
should  be  noted  that  discrepancies  in  amplitude  between 
direct  integration  and  the  Weber  function  techniques  outside 
the  interference  zone  may  be  overlooked,  since  certain 
approximations  used  in  developing  the  theory  required  the 
epicentral  distance  to  be  close  to  the  critical  point. 
Results,  displayed  in  Figure  3.6,  indicate  the  correctness 
of  the  Weber  function  approach  within  the  interference  zone. 


3 . 4  Conclusions 


Approximate  methods  of  determining  seismic  ray 
amplitude  have  recently  become  in  vogue  for  the  theoretical 
investigation  of  wave  propagation  in  homogeneous, 
pi ane- 1 ayered  media.  Two  of  these,  asymptotic  ray  theory  and 
the  Weber  function  approach  for  interference  reflected  head 
waves,  have  been  shown  to  match  very  closely  with  a  more 
exact  method  of  direct  numerical  evaluation  of  the 
displacement  potentials.  The  matching  between  the 
integration  method  and  the  two  high  frequency  approximations 
became  better  as  the  source  frequency  was  increased. 
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In  summary,  this  chapter  has  provided  a  numerical 
investigation  of  two  methods  of  seismic  ray  amplitude 
computation:  asymptotic  ray  theory  and  the  Weber  function 
approach.  These  methods  complement  each  other  in  reference 
to  their  regions  of  applicability.  This  renders  them 
especially  useful  for  purposes  of  constructing  synthetic 
sei smograms . 


4.  RAY  TRACING  BY  CLASSICAL  METHODS 


4 . 1  Introduction 


Classical  methods  of  computing  ray  paths  for  seismic 
energy  have  been  Known  for  some  time.  One  of  these  consists 
of  a  variational  approach  of  Fermat's  Principle.  We  shall 
term  this  technique  "The  Bending  Method".  Another  classical 
method  consists  of  numerically  solving  a  system  of  first 
order  differential  equations.  For  certain  simple  velocity 
distributions  (i.e.  simple  enough  to  be  specified  by  a 
certain  analytic  function  of  position)  ray  paths  may  be 
analytically  determined  by  complex  analysis.  The  remaining 
ray  tracing  methods  rely  upon  solving  a  set  of  coupled 
non-linear  differential  equations.  For  this  purpose, 
computationally  costly  numerical  methods  must  be 
implemented.  A  newer  method,  employing  local  velocity 
gradients  for  which  analytic  expressions  are  Known  will  be 
presented  in  the  next  chapter. 


4 . 2  The  bending  method  for  ray  tracing 

The  bending  method  as  described  by  Julian  and  Gubbins 
(1977)  was  considered  for  purposes  of  tracing  seismic  rays 
through  two-dimensional  media.  Additional  authors  include 
Wesson  (1970)  and  Chander  (1975).  The  method  basically 
perturbs  an  initial  guess  at  the  true  ray  path  between  the 
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source  and  receiver  until  a  stationary  time  path  is 
obtained.  The  differential  equations  for  the  ray,  obtained 
by  the  calculus  of  variations,  are  expressed  in  terms  of 
changes  in  the  ray  path  coordinates.  The  system  of  equations 
is  linearized  and  solved  by  standard  techniques.  This 
solution  represents  corrections  to  the  ray  path.  The  process 
is  reiterated  until  either  the  RMS  path  perturbation  between 
consecutive  iterations  is  less  than  some  prescribed  error 
bound  or  the  variational  equations  for  the  ray  (expressed  by 
finite  differences)  are  solved.  This  method  is 
computationally  expedient  if  an  efficient  method  is  used  to 
solve  the  linear  algebraic  equations.  For  this  purpose,  one 
may  take  advantage  of  the  fact  that  the  coefficient  matrix 
involved  is  tridiagonal  and  apply  Gaussian  elimination 
algorithms  as  ably  discussed  in  Henrici  (1962). 

Figure  4.1  gives  an  example  of  the  use  of  this 
technique.  The  solid  path  is  the  true  ray  path  for  the 
velocity  structure  v(z)=3+1.3z.  This  ray  is  an  arc  of  a 
circle,  the  equation  for  which  is  given  in  many  standard 
texts  (Nettleton  1940).  The  initial  guess  at  the  ray  path 
was  taken  to  be  a  straight  line  joining  source  and  receiver. 
The  corrected  path  is  shown  after  N=2,4,9  iterations.  We 
have  also  displayed  the  error  in  the  ray. 

For  purposes  of  ray  tracing  alone,  this  technique 
appears  to  be  of  great  service.  The  method  has  the  desirable 
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feature  of  allowing  the  user  to  specify  the  receiver  point, 
which  other  methods  seek  to  reach  by  varying  the  take-off 
angle  of  the  ray  at  the  source.  However,  for  purposes  of 
synthetic  seismogram  computation  at  several  receiver 
locations  in  a  complicated  model,  the  method  has  some 
drawbacks.  Let  us  define  a  family  of  rays  as  all  those  rays 
which  travel  outward  from  the  source,  strike  the  same 
interfaces  in  the  same  order ,  are  converted  from  P  to  S  mode 
or  vice  versa,  and  finally  arrive  at  the  free  surface.  Due 
to  a  complex  two-dimensional  velocity  structure  (Psencik 
1972)  or  curvilinear  interfaces  (Hron  1973)  or  both,  this 
family  of  rays  may  produce  a  large  number  of  individual 
travel  time  branches.  For  a  receiver  at  some  distance  from 
the  source  the  bending  method  would  have  to  be  utilized  to 
find  the  ray  path  corresponding  to  each  arrival  (i.e.  for 
each  branch  which  exists  at  the  distance).  This  would 
require  some  initial  guess  for  each  of  these  ray  paths.  It 
is  at  this  point  that  the  method  is  insufficient.  Firstly, 
the  ray  paths  in  complex  models  may  be  rather  exotic  (Hron 
et  al  1977)  and  an  initial  estimate  may  be  very  difficult. 
Secondly,  a  poorly  conjectured  initial  path  will  converge  to 
a  ray  lying  on  a  different  branch  than  the  one  desired. 
Finally,  this  method  does  not  give  information  on  the  total 
number  of  travel  time  branches  related  to  the  same  family  of 
rays,  nor,  more  crucially,  whether  a  particular  arrival  lies 
on  a  forward  or  reverse  branch.  This  information  is 
indi spensible  in  determining  the  phase  shift 
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tt/2  associated  with  the  geometrical  spreading  (Cerveny 
and  Ravindra  1971,  p.  76,  Choi  1978)  for  rays  having  passed 
through  a  number  of  caustics.  For  these  reasons,  we  must 
seek  alternative  approaches  for  seismic  ray  tracing. 


4 . 3  Complex  analysis 

A  rather  limited  approach  to  the  problem  was  presented 
by  Barnes  and  Solomon  (1971).  Their  technique  requires  an 
analytic  expression  for  the  velocity  field  of  the  medium. 
For  certain  fields,  ray  paths  can  be  derived  using  the 
theory  of  complex  variables.  However,  it  is  difficult  to 
find  a  function  which  would  even  partially  correspond  to  the 
velocity  distribution  of  a  real  medium.  For  this  purpose  we 
will  not  further  investigate  this  method. 


4 . 4  Pi f f erent i a  1  equations  for  ray  tracing 

Ray  paths  for  seismic  energy  may  be  computed  by  solving 
a  system  of  first  order  differential  equations.  Various 
authors  have  used  different  approaches  to  arrive  at  these 
equations,  including  Eliseevnin  (1964),  Chen  and  Ludwig 
(1972),  Psencik  (1972),  Green  (1976)  and  White  (1978),  among 
others . 


In  two-dimensional  media  with  a  point  source,  three 
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differential  equations  must  be  solved  for  the  coordinates  of 
the  endpoint  of  the  ray,  (x(t),  z(t)),  and  the  angle  that 
the  ray  makes  with  the  vertical,  0(t)  ,  as  a  function  of 
time,  t.  Following  Psencik  (1972)  these  may  be  written  as 

x  v  sinQ 

i  | 

i 
' 

d 

dt  z  i  =  |  v  cose 

9  i  v  sin0  -  v  cos6 

.J  i_  z  x 

where  v(x,z)  is  the  velocity  of  propagation  ( P  or  S  wave 

mode)  at  the  point  under  consideration,  v  ,  v  are  its 

x  z 

spatial  derivatives.  These  equations  must  be  solved  subject 
to  the  initial  conditions 

x ( 0 )  =  xQ,  z ( 0 )  =  Zq,  e(O)  =  eQ 

where  (xQ,zQ)  is  the  source  location  and  0Q  is  the  take 
off  angle.  Note  that  for  equation  4.1  to  hold  8  must  be 
measured  counter-clockwise  from  the  negative  z  axis  to  the 
tangent  of  the  ray  path.  We  take  the  positive  z  axis  to 
point  into  the  earth  from  the  free  surface. 

The  solution  to  equation  4.1  requires  knowledge  of  the 
velocity  field  and  its  spatial  derivatives  everywhere  that 
the  ray  travels.  The  input  velocities  may  be  specified  on  a 
rectangular  grid  rather  than  an  analytical  expression.  This 
is  due  to  the  difficulty  in  finding  a  function  which 
reasonably  approximates  true  velocity-position  relationships 


0  <6  <  2tt 


(4.1  ) 
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in  the  earth.  In  order  to  determine  v  and  its  derivatives 
within  any  rectangle  of  the  grid  one  may  use  linear 
interpolation  among  the  four  corner  points.  This 
approximation  was  utilized  by  Psencik  (1972).  However,  this 
will  lead  to  second  order  discontinuities  in  the  ray  from 
rectangle  to  rectangle.  An  improved  technique  has  been 
implemented  to  employ  bicubic  splines  for  the  velocity 
field.  By  this  method  the  velocity  at  any  point  is  dependent 
upon  not  only  the  velocity  at  the  four  nearest  grid  points, 
although  these  are  weighted  more  heavily,  but  also  by  the 
general  trend  of  the  velocity  in  the  entire  layer. 

Consider  a  point  (x,z)  which  lies  within  a  grid 
rectangle  with  (xi,zz  as  coordinates  of  its  upper  left 
corner.  By  employing  all  the  velocity-position  data  given 
within  the  layer  that  houses  (x,z)  we  may  compute  bicubic 
spline  coefficients  in  the  form  of  4x4  matrices,  C.^  , 
i ,j=1 ,2,3,4,  (Ahlberg  et  al  1967,  Chapter  7).  A  spline 
interpolation  may  then  be  made  as 

v(x,z)  =  Cij(x-xz)i  1(z-zz)J  1 

where  we  sum  over  i  and  j.  If  one  wishes  to  use  the  same 
velocity  model  for  many  rays  in  different  computer  runs,  a 
considerable  saving  of  central  processing  unit  (CPU)  time 
may  be  achieved  by  computing  all  spline  coefficients  once, 
writing  them  on  a  permanent  direct  access  file  in  the 
computer  memory,  and  subsequently  scanning  the  file  for  the 
appropriate  coefficients  as  the  program  progresses. 
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The  effect  of  geological  curvilinear  interfaces  where 
the  ray  directions  change  di scont i nuous ly  must  be  taken  into 
account  in  solving  equation  4.1.  By  the  principle  of  the 
isolated  element  and  Snell's  Law,  we  may  determine  the 
initial  parameters  for  solving  equation  4.1  after  the  ray 
impinges  an  interface.  The  location  of  interfaces  may  be 
spcified  by  a  series  of  (x,z(x))  coordinates.  We  may  smooth 
this  initial  data  by  the  use  of  cubic  splines  (Ahlberg  et  al 
1967,  Chapter  2).  This  approach  will  also  allow  us  to 
determine  the  interface  normals  quite  easily. 

Although  the  formulation  for  ray  tracing  presented  here 
will  give  us  ray  paths  which  can  be  considered  to  be 
correct,  it  has  a  serious  impediment.  This  is  the  extremely 
large  amount  of  CPU  time  needed  to  solve  the  differential 
equations  with  a  reasonable  degree  of  accuracy.  The  severity 
of  this  is  compounded  when  many  rays  are  considered,  as  in 
the  construction  of  synthetic  seismograms.  To  alleviate  this 
the  next  chapter  presents  a  ray  tracing  concept  which 
consumes  only  a  fraction  of  the  CPU  time  required  by  the 
differential  equation  process  and  still  yields  reasonably 
accurate  results. 
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5.  RAY  TRACING  BY  CIRCULAR  APPROXIMATION 


5 . 1  Introduction 


Chapter  4  described  techniques  for  tracing  rays  through 
two-dimensional  media  with  curvilinear  boundaries.  The 
differential  equation  method  is  precise  in  that  it  takes 
into  account  all  the  known  values  of  velocity  within  a  given 
layer  to  compute  the  velocity  at  any  point.  However,  from  a 
computational  viewpoint,  the  method  is  extremely  expensive. 
An  approximate  method  has  been  developed  which  gives 
reasonably  close  results  to  the  exact  theory.  Much 
computational  effort  may  be  spared  if  one  approximates  the 
velocity  field  v(x,z)  within  a  small  region  of  the  model  to 
be  expressible  as 

v(x,z)  =  v 0 (x q / z q )  +  grad  v-dr 

where  vQ  and  (grad  v)=K  are  constants.  With  this 
representation,  one  may  solve  classical  equations  to  direct 
the  ray  endpoint  to  its  proper  location  as  a  function  of 
time.  Similar  techniques  were  developed  independently  by 
Gebrande  (1976),  Will  (1976),  Aric  and  Gutdeutsch  (1978)  and 
Dantz  (1978).  All  of  these  techniques  require  the 
approximation  of  the  isovelocity  lines  of  the  medium. 

Rather  than  using  the  computationally  expensive  cubic 
splines  to  represent  a  geological  interface  from  a  set  of 
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discrete  data,  we  shall  approximate  the  interface  by  a  set 
of  piecewise  continuous  parabolae.  The  use  of  these  second 
degree  equations  will  allow  us  to  compute  analytical  rather 
than  numerical  solutions  for  the  intersection  between  a  ray 
and  an  interface. 


5 . 2  Ci rcular  approximation  ray  tracing  wi thout  interfaces 


Consider  an  elastic  isotropic  medium  whose  velocity 
field  v(x,z)  is  specified  by  velocities  given  on  the  nodal 
points  of  a  rectangular  grid.  The  spacing  between  grid 
points  may  vary  over  the  whole  grid.  Let  the  coordinate 
system  by  which  these  grid  points  are  given  be  termed  the 
home  coordinates.  Place  a  point  source  of  harmonic  waves 
somewhere  within  this  medium,  say  at  (  xs'zs  )■  Let  us 
define  another  set  of  coordinate  axes  parallel  to  the  home 
coordinates  and  with  origin  at  the  source.  This  system  of 
coordinates  will  be  termed  the  translated  coordinates. 


The  rectangle  in  which  the  source  lies  may  be  divided 
into  two  equal  triangles  with  the  source  lying  in  either  the 
upper  or  lower  triangle.  This  decomposition  of  the  model' s 
blocks  into  triangles  is  the  basic  advantage  of  this  ray 
tracer,  for  we  may  employ  velocities  at  the  three  corner 
points  of  the  triangle 
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vi  (x± '  zj_  ^  ,  i  =  1 ,2,3  to  compute  a  unique  linear  mapping 


v(x,z)  =  Ax  +  Bz  +  C 


(5.1  ) 


by  solving  three  linear  algebraic  equations  for  the  three 
unknown  coefficients  A,BtC.  Numerical  experimentation 
demonstrated  that  it  is  most  beneficial  to  divide  the 
rectangle  by  the  diagonal  joining  the  two  points  adjacent  to 
the  corner  point  of  the  rectangle  nearest  to  the  source 
point . 

Figure  5.1  shows  a  typical  geometry  for  this  situation. 
If  the  local  velocity  field  is  given  by  equation  5.1,  we  do 
not  really  have  a  depth-dependent  velocity  as  pledged. 
However,  if  the  coordinate  axes  are  rotated  through  an  angle 


(5.2) 


then  we  can  define  a  set  of  translated,  rotated  coordinates 
(x 1 , z ' )  such  that 


!  x 1  cos,  sincj>  x 


z 


-sing  coscj)  z 


(5.3) 


In  Figure  5.1,  note  that  is  positive  when  measured  in 

the  clockwise  direction.  In  the 
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Figure  5.1  Symbols  and  geometry  used  to  compute  raypaths  by 
circular  approximation. 
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(x',z')  coordinates  the  velocity  field  is  given  by 


vU\z')  =  v(0,0)  +  (A2+B2^z' 


(5.4) 


Nettleton  (1940)  has  derived  analytic  expressions  for 
seismic  ray  paths  for  a  medium  whose  velocity  field  is  given 
by  v  =  vQ  +  kz .  Essentially,  these  paths  are  circles  in  the 
(x',z!)  plane  defined  by 
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ksmfi 


0 


(5,5) 


where  0O  is  the  take  off  angle  (measured  from  the 
vertical)  of  the  ray,  0<9  <2tt  .  As  well,  the  travel  time  T 
to  reach  any  point  (x',z’)  on  the  ray  path  is 


1  ,  -1 
T  =  t  cosh 
k 


k2  ( x  * 


2  ,  2 . 
+  z '  ) 


2vv 


+  1 


0 


(5.6 


Of  course,  ifv1  =  v9  =  v3,  the  block  is  homogeneous  and  the 
rays  degenerate  to  straight  lines. 

The  ray  can  now  be  traced  through  the  triangle  up  to 
its  exit  point.  Determination  of  the  exit  point  can  be  posed 
as  the  problem  of  finding  the  intersection  between  a  circle 
(equation  5.5)  and  a  straight  line 
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z'  =  Mx'  +  B.  These  straight  lines  are  the  lines  bounding 
the  triangle  in  the  translated  rotated  coordinates.  We  must 
find  the  point  of  intersection  ^xe'ze^  and  the 
correspondi ng  time,  .  Since  there  may  be  up  to  six 

intersection  points,  the  criteria  for  determining  the 
correct  one  are: 

(1)  those  for  which  (x',z')  lies  within  the 

0 

specified  endpoints  of  the  lines  z  =  Mx '  +  b 

(2)  that  for  which  the  time  ,  is  minimal. 

The  ray  may  either  have  entered  the  second  triangle  of 
this  block  (by  crossing  line  2)  or  a  completely  new  block 
(by  crossing  lines  1  or  3).  If  the  former  event  results,  the 
procedure  is  repeated  starting  with  the  calculation  of 
coefficients  A,B,C.  If  the  latter  presides,  the  operation  is 
initiated  beginning  with  division  of  the  new  block. 

This  technique  of  ray  tracing  requires  considerably 
less  CPU  time  than  the  exact  differential  equation  method 
described  in  Chapter  4.  The  chief  saviour  is  the  use  of 
analytic  expressions.  However,  one  difficulty  must  be 
highlighted.  Since  the  velocity  field  is  only  specified 
locally,  the  ray  will  "see"  the  lines,  separating  the  blocks 
as  second  order  interfaces  resulting  in  a  discontinuous 
derivative  of  the  ray  angle  upon  crossing  these  lines. 
Numerically,  the  effect  is  minimal,  as  the  next  chapter  will 
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bear  out. 


5 . 3  Determining  the  ray  as  a  function  of  time 


For  purposes  of  producing  ray  diagrams  it  may  be 
beneficial  to  Know  the  location  of  the  ray  endpoint  and  the 
angle  with  respect  to  the  rotated  vertical  as  a  function  of 
time.  The  time,  t,  is  measured  from  the  moment  the  ray 
enters  the  triangle.  To  disclose  these  quantities  one  must 
solve  equation  4.1  with  a  velocityv  =  vQ  +  kz  and  x(0)  =  x 


0  ' 


z (0)  =  Zq'  u vw/  ~o 


0(0)  =  @n  .  The  solutions  to  these  equations  are 


v, 


0 


x(t)  =  tan  (e2kt-l)  (l+ae2kt)  X 


z(t)  =  ^  (ekt(l+a) (l+ae2kt)_1  -1) 


_  X  kt  0 

0(t)  =  2tan  (e  tan  -~-) 


(5.7) 


where 


v(t)  =  v0(l+a)ekt(l+ae2kt)  1 


a  =  tan 


2  90 


5.4  The  effect  of  an  interface 


The  approximation  of  interfaces  by  piecewise  continuous 
parabolae  is  discussed  in  Appendix  C.  If  curvilinear 
interfaces  are  to  be  introduced  into  the  model,  record  must 
be  Kept  of  which  blocKs  are  cut  by  these  interfaces.  If  the 
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ray  enters  a  block  which  contains  an  interface,  specific 
steps  must  be  taken.  The  first  of  these  is  to  eliminate  the 
velocity  datum  on  any  corner  point  of  the  block  which  lies 
on  the  other  side  of  the  interface  from  which  the  ray  is 
incident.  The  datum  must  be  replaced  by  extrapolation  among 
the  velocities  at  the  points  nearest  the  point  in  question. 
In  practice,  this  can  be  performed  by  linear  extrapolation 
utilizing  the  three  nearest  points. 

The  second  modification  required  in  blocks  cut  by  an 
interface,  is  to  determine  the  location  and  time  of  the 
ray's  intersection  with  the  interface  (if  any).  If  we  have  a 
parabola  described  as  z ( x ) = ( ax+b ) x+c  and  a  circle  given  by 
( x-A ) ( x-A ) + (z-B ) (z-B ) =r*r  then  the  (at  most  4)  points  of 
intersection  are  the  roots  of  the  following  fourth  degree 
polynomial 

a2x^ +  2abx2  +  (b2+2a (c-B) +1) x2 

+2 (b (c-B) -A) x  +  (A2-r2 +  (c-B)2)  =  0 

The  equations  for  the  circle,  parabola  and  intersection 
points  are  to  be  taken  in  the  unrotated,  translated 
coordinates  with  the  origin  at  the  point  where  the  ray 
enters  the  triangle  containing  the  interface.  The  fourth 
degree  polynomial  can  be  solved  numerically  or  sometimes 
even  analytically.  The  requirements  for  an  intersection 
point  to  be  the  true  one  are  the  same  as  for  the 


■ 


88 


ray-triangle  boundary  case. 

Using  parabolae  to  approximate  interfaces  makes  i t 
simple  to  compute  the  interface  normals  at  the  point  of 
incidence.  Knowing  these,  as  well  as  the  incident  angles  at 
the  time  of  intersection  (equation  5.7),  it  is  an  elementary 
operation  to  reflect  or  refract  rays  according  to  Snell's 
Law. 


The  next  chapter  will  discuss  ways  that  this  method,  as 
well  as  the  method  of  differential  equations,  can  be 
utilized  to  calculate  ray  amplitude  at  the  earth's  surface. 


6.  GEOMETRICAL  SPREADING  OF  WAVEFRONTS  IN  LATERALLY 


INHOMOGENEOUS  MEDIA 


6 . 1  Introduction 

Computation  of  ray  paths  in  laterally  inhomogeneous 
media  has  been  discussed  in  the  preceeding  chapters. 
Applicable  methods  for  this  task  include  solution  of  a 
system  of  differential  equations  and  the  method  of  circular 
approximation.  The  problem  of  determination  of  geometrical 
spreading,  and  hence  amplitude,  is  solvable  by  techniques 
similar  to  those  for  ray  paths.  In  this  chapter  equations 
will  be  given  and  numerical  tests  performed  to  obtain  ray 
amplitudes  and  travel  times. 


6.2 


D i f ferent i a  1 


equation  solution  for  ray  amp! i tude 


Differential  equations  have  been  derived  by  Cerveny  et 
al  (1974),  among  others,  by  which  the  cross  sectional  area 
of  the  ray  tube  may  be  expressed.  A  brief  summary  of  these 
equations  will  be  given  here.  For  more  in-depth  discussion, 
the  reader  is  referred  to  the  original  publication. 
Basically,  the  equations  are  an  extension  of  the  ray  tracing 
equations  discussed  in  Chapter  4. 


We  introduce  ray  coordinates 


. 
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ql'g2:>T  where  <5i’g2  specify  the  ray  under 
consideration  and  t  the  position  of  a  point  along  the  ray. 
Physically,  ^  may  represent  the  declination  of  the  ray  at 
take-off  and  the  azimuthal  angle.  Thus  any  point  along 
the  ray  may  be  described  by  x  =  (x^ (q1,q2 , t) )  f  i  =  1  ,2,3. 
The  cross  sectional  area  of  an  elementary  ray  tube,  dcJ  can 
be  evaluated  by 

da  =  J  dq1dq2  ,  > 


where 


5x  x  3x 
3q1  9q2 


(6.2) 


and  x  =  (x^x^x^  .  In  zeroth  order  ray  theory,  the 

amplitude,  A,  of  a  wave  at  time  T  is  expressible  as 


A(t1)  = 


v(To>  ety  J(Toi  p 

vfTj^)  p  (  T  ^  )  J  (  T  ) 


A(To} 


(6.3) 


where  v  denotes  velocity  and  p  density.  T0  is  some 
reference  time.  It  is  convenient  to  take  this  time  to 
correspond  to  the  wave  front  on  the  unit  sphere  about  the 
source.  Within  this  sphere,  the  medium  is  assumed 
homogeneous.  The  components  of  the  slowness  vector 
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(pi(T))  can  be  written  in  terms  of  the  azimuth  angle  A 
and  dec  1 i nat ion  D 

_  cosAsinD  ...  sinAsinD  cosD 

P1 - v - '  P2 - 7 - '  P3  =  -T- 


Let  us  introduce  the  following  notation: 


yai(ql' 


q2'T)  = 


9xi{q1,q2/T) 


3q 


a 


z  . 
ai 


t) 


3pi (q1/q2f T) 


3q 


a 


(6.4) 


(6.5) 


for  a  =  1,2  .Using  the  Laplace  identity,  namely 

(a*b) • (c*d)  =  (a*c) (b*d)  -  (b*c) (a-a) 

we  can  rewrite  6.2  in  a  more  useful  form 


J  = 


(6.6) 


where 

E  =  yliyli'  F  =  y2iy2i'  G  =  yliy2i 

and  summation  convention  is  assumed.  Functions  are 

obtained  by  means  of  partial  differentiation  of  equation  4.1 
wi th  respect  to  q^  .  These  are 


dy 


oti 


dx 


p .  y  . 

^1  CLJ 


2v 


3v 
3x  . 


+  v 


ai 


^2ai  _  j  3v  3v 

2  3x .  3x . 
v  _  i  3 


3  2v 


3x . 3x . 
i  3  _ 


dx 


v 


(6.7) 
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Similar  equations  to  these  are  available  in  Wesson  (1970)  or 
in  Green  (1976)  for  spherical  geometry.  The  system  must  be 
solved  simultaneously  with  the  ray  tracing  equations,  4.1, 
to  obtain  both  Kinematic  and  dynamic  properties  of  waves  in 
inhomogeneous  media.  The  numerical  problems  in  solving  this 
system  are  the  same  as  those  discussed  in  Chapter  4.  The 
user  can  now  expect  a  substantial  increase  in  CPU  time  to 
solve  the  extra  equations. 


6 . 3  The  effect  of  an  interface 

The  differential  equations  4.1  and  6.7  can  be  solved  in 
a  region  where  v(x^)  is  a  smooth  function  of  the 
coordinates,  i.e.  up  to  the  point  where  the  ray  intersects 
with  an  interface  of  first  order .  We  assume  that  this 
interface  is  given  by  an  equation  f(x^)  =  0  ,  where  f  is  to 
have  continuous  partial  derivatives  up  to  second  order  in 
some  neighborhood  of  the  intersection  point.  The  components 
of  the  unit  normal  are 

N.  =  e  liL/r^^-T5 
1  dx.  /  9  X  ■  9x  . 

1  _  D  HJ 

where  e=±l  ,  chosen  so  that  the  normal  vector  points  into 
the  layer  from  which  the  ray  is  incident.  At  such  interfaces 
the  quant i ties 
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change  di scont i nuous ly .  Cerveny  et  al 


(1974)  have  presented  equations  for  the  change  in  ^ai'zai 
while  the  change  in  is  obtainable  through  Snell's 

Law.  Let  quantities  which  are  to  be  evaluated  at  the 

interface  be  written  in  upper  case.  Variables  with  a  tilde 

represent  the  value  on  the  reflected  or  refracted  side  of 

the  interface  while  non-tilded  variables  are  measured  on  the 
incident  side.  The  reflected  side  coincides  with  the 
incident  side  of  the  interface,  in  that  energy  is  carried 
away  from  the  interface  into  the  same  medium  from  which  it 
arrived.  According  to  equations  38  and  39  of  the  above  paper 


(6.8) 


Nk  3V  3T  _  X  V3V 


V3V 


a 


(6.9) 


where 
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0  <A  =  P.N. 

l  l 

B  =  (V“2  -  V"2  +  A2f*>0 


The  upper  sign  is  used  for  a  refracted  wave  and  the  lower 
for  a  reflected  one.  The  derivative  of  an  arbitrary  function 
with  respect  to  along  the  interface  is 


3G 

3q 


a 


3G 

3x 


Y 


l 


ai 


+ 


v2p 


3  T 


i  3q 


a 


To  determine  the  amplitude  at  a  time  along  a  ray 
that  has  encountered  K  interfaces  we  employ  the  following 
equation  from  asymptotic  ray  theory 


rv(V  p(V  J(V 


A(T1]K=  A(T0)  I  J(TX^ 


k 

n 

j=i 


v(t.)  p  (t  .  )  J  (t  .  )  % 

3 _ 2 —  - 2 —  r  . 

v(T.j)  P  ( T  j  )  J(Tj)  3 


(6.10) 


where  Rj  are  the  correspondi ng  plane  wave 
refraction  or  surface  conversion  coefficients 
and  Ravindra  1971)  and  are  the  times  of 

with  the  interfaces.  It  should  be 

J(Tj)/J(T^)=  cos  1/ cos  i  where  i  (  I  )  is 

incidence  (reflection  or  refraction)  Cerveny  et 


reflection , 

( see  Cerveny 
intersection 
noted  that 
the  angle  of 
al  (1977). 


Thus,  the  computation  of  rays  and  ray  amplitude  in 


laterally  inhomogeneous  media  can  be  performed  by  the 
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following  procedure: 

(1)  Integrate  equations  6.7  and  4.1  up  to  the 
point  of  intersection  with  an  interface.  In 
this  fashion  we  obtain  P.  ,Y  .  ,z  . 

i  ai  oti 

(2)  From  the  relations  6. 8-6. 9  obtain  values  for 

Yai'^cti-  is  obtained  from  Snell's  Law. 

(3)  At  the  interface  compute  the  product  term 
required  in  equation  6.10. 

(4)  Use  the  values  of  Pi'Yai'Zai  as  new  initial 
values  and  begin  solving  equations  6.7  and  4.1 
subject  to  the  new  initial  values. 

(5)  Steps  ( 1 ) - ( 4 )  are  repeated  until  either 
another  interface  is  reached,  the  ray  is 
terminated  at  a  receiver  or  the  ray  is  traced 
out  of  bounds. 

Numerical  examples  of  this  method  will  be  discussed 
further  in  this  chapter. 


6 . 4  Ray  amp  1 i tude  by  dynami c  ray  tracing 


In  a  recent  work,  Cerveny  and  Hron  (1979)  studied  the 
problem  of  computation  of  ray  amplitude  in  a  generally 
inhomogeneous  medium  with  curvilinear  interfaces.  They 
introduced  the  process  of  "dynamic  ray  tracing",  which 
employs  second  derivatives  of  the  seismic  time  field  with 
respect  to  the  ray  coordinates.  For  a  laterally 
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inhomogeneous  two-dimensional  medium,  the  dynamic  ray 
tracing  system  reduces  to  two  non-linear  differential 
equations.  The  phase  matching  method  is  used  to  determine 
discontinuities  of  individual  quantities  in  the  system  when 
the  wave  impinges  a  curved  boundary  separating  two  generally 
inhomogeneous  media.  The  final  equations  only  and  some 
numerical  examples  will  be  presented  in  this  chapter  of  the 
thes i s . 

In  a  two-dimensional  model,  in  which  elastic  and 
geometrical  characteristics  do  not  change  in  one  fixed 
direction,  the  propagation  of  seismic  waves  can  be 
investigated  in  the  plane  perpendi cu 1 ar  to  this  fixed 
direction.  Under  these  circumstances,  the  dynamic  ray 
tracing  system  of  Cerveny  and  Hron  (1979)  reduces  to 

dM  =  v2  _  M?  \j  i  \ 
dT  v 


subject  to  M ( 0 )  =  N ( 0 )  =  0 .  Equation  6.11  corresponds  to 
equation  93  of  the  above  reference.  The  variable  M 
represents  the  product  of  the  wavefront's  radius  of 
curvature  (in  the  plane  of  propagation)  multiplied  by  the 
velocity  of  propagation  at  that  time, 
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t  .  The  variable  N  is  similar  but  for  the  radius  of 
curvature  perpendicular  to  the  plane  of  propagation.  VII 
represents  the  second  derivative  of  velocity  in  the 
direction  tangent  to  the  wavefront. 


Unfortunately,  the  spreading  function,  J,  employed  in 
equation  6.10  for  ray  amplitude  is  not  readily  available 
from  M  and  N.  It  is  given  by 


J  (  T  ) 


T 


J(tq) 


exp 


2  1  ^  1  .r 

V  M  +  N  dS 


0 


(6.12) 


If  the  first  ray  segment  emanating  from  the  source  is  under 

consideration,  then  =  Vq1  ,  corresponding  to  a 

homogeneous  unit  sphere  surrounding  the  source.  vQ  is  the 

velocity  at  the  source  and  J(t  )  =  sin0  where  e  is  the 

o  0 

ray  take  off  angle  as  measured  from  the  vertical  at  the 
source.  Otherwise,  To  is  the  time  at  the  beginning  of  the 
current  ray  segment. 


When  the  ray  impinges  a  curved  interface  between  two 
different  media,  the  quantities  M  and  N  in  equation  6.11 
will  become  discontinuous.  The  initial  values  of  M  and  N  on 
the  reflected  or  refracted  side  of  the  interface,  denoted  by 
a  tilde,  must  be  computed  in  order  to  begin  solving  equation 
6.11  again.  According  to  the  theory  of  dynamic  ray  tracing, 


M 


.  2 
COS  1 

cos  i 


1 

M 


+  K  + 


I 


-1 


' 
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N 


=  N 


where 


K  = 


uC 


(cos  i) 


(6.13) 


I  = 


sin  i 


— 

\vcos 

\ 


9_v 

9x 


sin 


i 


cos 


i 


8v 

-  7T—  sin 

dX 


i 


7T—  cos  i 

d  Z 


cos  i  ,  cos  i 


i (i)  is 

the  positive 

acute 

angle 

that 

the 

i ncident 

( reflected 

or  refracted) 

ray 

makes 

wi  th 

the 

interface 

normal.  v(v)  is  the  velocity  of  the  incident  (reflected  or 
refracted)  wave  at  the  incident.  In  u,  the  upper  (lower) 
sign  is  taken  for  reflected  (refracted)  waves.  C  is  the 
curvature  of  the  curved  boundary,  which  is  taken  as  positive 
(negative)  for  those  parts  of  the  boundary  which  are  seen  as 
convex  (concave)  by  an  observer  located  in  the  medium  from 
which  the  ray  is  incident. 


K  is  termed  the  curvature  term,  describing  the  effect 
of  interface  geometry  upon  ray  amplitudes.  On  the  other 
hand,  I,  depending  on  variations  of  the  elastic  parameters 
near  the  interface,  is  termed  the  inhomogeneity  term.  This 
is  a  key  advantage  of  dynamic  ray  tracing,  in  that  it  allows 
one  to  study  separately  the  effect  of  inhomogeneity  and 
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interface  geometry  upon  ray  amplitude.  Numerical  examples  of 
this  will  be  presented  later  in  this  chapter. 


6 . 5  Ray  amp! i tude  by  ci rcu 1  a r  approximation 


The  amplitude  of  a  wave  propagating  in  an  arbitrarily 
inhomogeneous,  isotropic  medium  is  given  by  (Cerveny  and 
Ravindra  1971,  equation  2.107) 


A  (M) 


q(ql'q2>  r  V0P0 
L  j  V  (M)  p  (M) 


K 

n 

j=i 


v(o. )p(o.)  h 
vToTTpToTy  Rj 


(6.14) 


where 

g 


M 

Mo 

6  (0  ) 


directivity  function  of  the  source  (taken  as  unity 
for  this  study ) , 

point  of  incidence  of  the  ray  upon  the  j-th 

interface  that  it  encounters, 

wave  velocity  and  density  at  source, 

receiver  location, 

some  reference  location, 

angle  of  incidence  upon  the  j-th  interface, 


L 


da (M)  \ 

da (Mq ) 


K 

n 

j=i 


da  (CK  )  N \\ 
da (0. ) 


The  product  in  the  geometrical  spreading  function,  L,  can  be 
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evaluated  as  the  ratios  of  cosines  of  angles  of  incidence 
and  reflection  (or  refraction)  as  in  Cerveny  et  al  (1974). 
This  term  represents  a  compensation  for  the  di scont i nui ty  of 
geometrical  spreading  across  the  interfaces.  The  tilded 
quantities  are  to  be  evaluated  on  the  reflected  or  refracted 
side  of  the  interface.  All  quantities  in  equation  6.14  are 
easily  computable  except  da(M)  .  We  shall  develop  an 
algorithm  to  compute  da(M)  for  laterally  inhomogeneous 
media  with  curvilinear  interfaces  by  means  of  the  circular 
approximation  discussed  in  Chapter  5. 

From  the  geometry  of  Fiaure  6.1 

da  (M)  =  cose  (M)  drd  pQ 

=  lr0cose(M)  deo  d<t,o 

where  r  is  the  magnitude  of  the  vector  pointing  from  source 
M0  to  receiver  M.  Within  the  unit  sphere  surrounding  the 
source  we  assume  the  model  to  be  homogeneous  allowing  us  to 
express  da(MQ)  =  sineQ  d0Q  d<J>g  .  Thus 


da  (M)  _  5r  cos9  (M)  8y 

da(MQ)  86q  sineQ  34>Q 


Numerical ly , 
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surface  at  M 

Figure  6. 1  Geometry  of  the  ray  tube  at  source  and  receiver 
The  heavy  lines  represent  the  earth's  surface.  The  angle 
0 1 s  confined  to  the  X-Z  plane,  while  a  is  confined  to 
the  Y-Z  plane.  v0 
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dr/dbO  may  be  determined  at  point  M  by  the  following 
procedure : 

(1)  Shoot  three  rays  on  either  side  of  the  ray 

arriving  at  the  receiver.  This  may  be  done  by 

subtracting  or  adding  small  amounts  A6Q  (  £&0 
-4 

^10  radians)  from  the  take-off  angle.  The  change 
in  eo  should  yield  a  change  in  epi central 
distance . 

(2)  Fit  a  cubic  spline  through  the  seven  (6  fr) 
data . 

(3)  Evaluate  the  derivative  3r/90o  at  the  take-off 
angle  which  causes  the  ray  to  arrive  at  the 
receiver  M. 


According  to  Cerveny  et  al  (1974),  g  =  3 y/8 4> 0  is 
defined  by  the  following  differential  equation  for  the 
two-dimensional  medium  which  is  under  consideration  in  this 
study: 


dg  = 
dt 


sin6 


v 


0  2  ,4.x 

-  V  (t) 


0 


Furthermore,  g  is  invariant  at  any 
Integrating  this  and  summing  over 
ray  passes  through  we  have 


sin0 


9  = 


0 


v 


0 


N 

l 

i=l 


2  2 
V0i  (1+ai> 


(zero  strike)  interface, 
all  N  triangles  that  the 


t . 
i 


2ki(t-ti-i> 


ti-i 


2k^ -2 


1+a  .  e 

l 


dt 


- 


■ 
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where 


...  velocity  at  entry  point  of  i-th  triangle, 


=  tan' 


Oi 


2 

angle  of  incidence  (in  rotated  coordinates) 

entry  of  i-th  triangle, 

velocity  gradient  in  i-th  triangle, 

total  ray  travel  time  upon  entry  (exit)  to 

i-th  triangle. 


upon 


( from) 


Note  that  tQ  =  0  and  tN  =  total  travel  time  for  the  ray. 
Upon  performing  the  integration  and  inserting  the  result  in 
equation  6.15  we  have 


da  (m) 

da (MQ) 


3r  cosG (M) 


36 


0 


v 


0 


N 

I 

i=l 


(l+ai) 

~2k~. 

l 


e  -1 

\ 


1+a  .  e 

l 


2k .  (t . -t .  .  )  -1 

l  l  i-l 


(6.16) 


If,  in  any  of  the  triangles,  k^  =  0  the  term  to  be  added 

2 

must  be  V0i(ti  ti-l)  The  computation  of  ray 

amplitudes  in  two  dimensionally  inhomogeneous  media  with 
curvilinear  interfaces  is  now  given  by  equations  6.14  and 


6.16. 
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6 . 6  Numer ica 1  tests 


FORTRAN  programs  have  been  written  to  compute  Kinematic 
and  dynamic  characteristics  of  waves  in  two-dimensional 
inhomogeneous  media  by  solving  differential  equations  and  by 
the  circular  approximation.  In  this  section,  numerical 
comparisons  will  be  presented  between  these  methods.  In  all 
these  models,  the  shear  wave  velocity  amd  volume  density  are 
given  by 

Vg(x,z)  =Vp(x/z)/3^p(Xf2)  =  3^  Vp  (x,z) 

from  Gardner  et  al  (1974). 

As  a  first  model  we  present  a  two-layered  laterally 
inhomogeneous  medium  with  curved  interfaces.  Since  this  is  a 
mathematically  contrived  model,  the  P  wave  velocity-position 
function,  V(x, y) ,  and  the  interfaces  are  given  analytically. 
Figure  6.2  is  a  description  of  the  model  showing  P1P2P2P1 
ray  paths  and  amplitude  distance  characteristics.  These 
amplitudes  were  computed  by  three  different  means,  namely, 
the  differential  equations  of  Cerveny  et  al  (1974),  the 
dynamic  ray  tracing  (D.R.T.)  system  of  Cerveny  and  Hron 
(1979)  and  the  circular  approximation  discussed  in  this 
chapter.  Results  are  given  as  a  per  cent  of  curve  1.  The 
differential  equations  technique  and  dynamic  ray  tracing 
both  required  about  two  minutes  of  CPU  time,  while  the 
circular  approximation  only  consumed  about  20  seconds.  With 


- 


105 


1  AMPLITUDE  BY  DIFFERENTIAL  EQUATIONS  ( CERVENY  .  LANDER.  PSENCIK) 

2  AMPLITUDE  BY  DYNAMIC  RAY  TRACING  (CERVENY  AND  HRON )  . . 

3  AMPLITUDE  3Y  CIRCULAR  APPROXIMATION  ( MARKS ) •— * . . 


110 


Pinnre  6.2  P1P2P2P1  in  a  mathematically  contrived  medium 
The  amplitudes  are  given  as  a  per  cent  of  curve  1. 


(INCLUDING  INTERFACE  CURVATURE  AND  INHOMOGENEITY  NEAR  INTERFACES  > 
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3 


amplitude  by 

AMPLITUDE  BY 
AMPLITUDE  BY 


DIFFERENTIAL  EQUATIONS  :  CERVENY  .  LANDER.  PSENC I K  )-«-*- 

DYNAMIC  RAY  TRACING  t CERVENY  AND  HftDN I - 

CIRCULAR  APPROXIMATION  I  MARKS  I 


Figure  6.4  Rays  P1S2S2P1  propagating  in  a  plane  layered 
medium.  The  seismic  wave  velocities  are  given  by  the 
mathematical  equations  presented  in  Figure  6.16. 
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this  saving  in  computational  time,  the  few  per  cent 
discrepancy  is  quite  acceptable. 

Dynamic  ray  tracing  allows  one  to  study  the  effect  upon 
amplitude  of  interface  curvature  and  velocity  inhomogeneity 
near  the  interfaces.  This  is  achieved  by  setting  either  K=0 
or  1=0  in  equation  6.13.  Figure  6.3  displays  the  result  of 
this  experimentation  for  the  model  shown  in  Figure  6.2. 
Setting  K=0  is  equivalent  to  considering  a  plane  interface 
at  the  point  of  incidence.  This  leads  to  a  15-22%  change  in 
amplitude.  Setting  1=0  is  equivalent  to  considering  a 
homogeneous  medium  surrounding  the  interface  at  the  point  of 
incidence.  This  yields  only  a  5%  change  in  amplitude. 
Clearly,  for  the  model  shown,  interface  curvature  plays  a 
Key  role  in  determining  seismic  ray  amplitude.  To  examine 
this  phenomenon  more  closely,  the  two  curved  interfaces  can 
be  replaced  by  flat  interfaces  at  depths  of  12  and  36  Km  as 
in  Figure  6.4.  To  demonstrate  the  programs'  capability  to 
handle  converted  phases,  results  are  displayed  for  the  ray 
P1S2S2P1.  With  the  removal  of  curved  interfaces,  the 
discrepancy  between  the  results  computed  by  three  different 
methods  decreases. 

The  Bohemian  Massif  is  a  laterally  inhomogeneous  medium 
taKen  from  a  seismic  survey  of  profile  VI  in  Czechos lovaKi a 
(PsenciK  1972).  Figure  6.5  exhibits  the  P  wave  isovelocity 
contours  for  this  model.  The  lateral  inhomogeneity  of  the 
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model  is  well  illustated  in  the  transect  plot  of  the 
velocity  field  shown  in  Figure  6.6. 

Figure  6.7  illustrates  a  comparison  between  ray  paths 
computed  for  identical  take  off  angles  at  the  source  (155.2 
Km) .  The  dashed  ray  segments  were  calculated  by  the 
differential  equation  approach,  while  the  solid  ones  were 
computed  by  the  circular  approximation.  Table  6.1  is  a 
detailed  account  of  the  results  of  this  test.  Note  that  an 
additional  11  seconds  of  computer  time  were  required  to 
evaluate  the  bicubic  spline  coefficients  for  the  velocity 
field.  This  was  performed  on  a  prior  computation. 

Figure  6.8  shows  a  comparison  between  the  travel  times 
for  the  ray  PI  PI.  In  Figure  6.9  we  have  compared  the 
logarithm  of  the  vertical  amplitude  (excluding  reflection 
coefficient).  Here  we  see  some  discrepancy  due  to  the 
different  methods  of  computing  geometrical  spreading,  L. 
However,  since  the  deviation  is  roughly  of  the  same 
magnitude  as  errors  in  determining  the  velocity  from 
inversion  of  seismic  data,  the  technique  of  circular 
approximation  appears  to  be  a  viable  instrument. 
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Comparison  of  P  Rays  Computed  by  Two  Methods 

Source  at  155.2  km 
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EP  I  CENTRAL  DISTRNCE  (KM) 

Figure  6.9  Comparison  of  amplitudes  (excluding  reflection 
coefficient)  for  the  ray  P1P1  as  computed  by  the  methods  of 
differential  equations  and  circular  approximation. 
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7.  NUMERICAL  MODELLING  OF  SEISMIC  BODY  WAVES  IN  LATERALLY 

INHOMOGENEOUS  MEDIA 


7 . 1  Introduction 


Methods  of  tracing  seismic  rays  and  computing 
geometrical  spreading  of  wave  fronts  in  two-dimensional 
laterally  inhomogeneous  media  have  been  discussed  in 
previous  chapters.  In  this  chapter,  a  general  procedure 
linking  all  these  processes  together  to  produce  synthetic 
seismograms  will  be  presented.  At  the  hub  of  this  package  of 
computer  programs  is  the  rapid,  yet  accurate  circular 
approximation  method  of  seismic  ray  tracing. 

7 . 2  The  relationship  between  take  off  anq 1 e ,  epi central 
di stance ,  travel  t ime  and  amp  1 i tude  i n  a  comp  1  ex  med i urn 

It  is  well  known  that  in  a  complex  geological  medium, 
the  travel  time  curve  of  any  phase  may  be  of  a  very 
complicated  nature.  Inhomogeneities  or  curved  interfaces 
lead  to  discontinuities  or  cusps  in  the  travel  time 
branches.  Triplication  or  even  higher  order  multiplicity  of 
travel  time  is  possible.  If  one  hopes  to  compute  ray 
theoretical  seismograms  then  a  necessary  requirement  is  to 
determine  the  ray  parameter  and  arrival  time  at  the  receiver 
of  all  phases  of  any  ray.  One  may  then  iterate  this  process 
with  as  many  rays  as  necessary  to  provide  a  reasonable 
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approximation  to  the  total  wave  field. 

Special  care  must  be  taken  to  locate  the  exact  position 
of  cusps  and  discontinuities  of  individual  travel  time 
branches  related  to  the  same  ray.  Since  we  are  dealing  with 
flat  Earth  models  of  finite  horizontal  extent,  it  is 
necessary  to  specify  left  and  right  boundaries  of  our  model. 
Thus,  rays  will  be  traced  out  of  the  bounds  of  the  model  for 
some  take  off  angles  at  the  source.  Let  us  introduce  T1  and 
T2  as  two  take  off  angles  which  cause  the  rays  associated 
with  them  to  be  traced  out  of  the  bounds  of  the  model. 
Furthermore,  assume  that  all  take  off  angles  between  T1  and 
T2  are  associated  with  rays  which  emerge  at  the  Earth's 
surface  within  the  bounds  of  the  model.  Then  the  travel 
time-distance  relationship  formed  for  angles  between  T1  and 
T2  is  defined  as  a  branch  of  the  total  travel  time  curve  for 
this  ray.  A  branch  may  consist  of  several  branch  parts  in 
which  epicentral  distance  is  a  monotonic  function  of  take 
off  angle  at  the  source.  Physically,  the  point  of 
intersection  between  two  adjacent  branch  parts  may 
correspond  to  either  a  caustic  or  a  grazing  ray.  An 
intelligent  ray  tracing  technique  should  be  able  to 
distinguish  which  is  correct.  As  a  simple  schematic  of  this 
situation  consider  Figure  7.1.  The  portion  of  the  curves 
from  A  to  B  corresponds  to  reflection  from  an  interface, 
while  segments  B  to  F  are  related  to  the  refracted  phases  of 
the  same  ray.  For  some  epicentral  distances  up  to  four 
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Figure  7.1  Take  off  angle,  epicentral 
and  amplitude  in  a  complex  medium. 
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arrivals  will  contribute  to  a  synthetic  seismograms.  For  a 
complex  geological  model,  the  number  of  individual  travel 
time  branches  as  well  as  the  number  of  monotonic  parts 
contained  within  each  branch  may  rise  to  rather  large 
values.  In  order  to  divide  the  individual  branches  into 
their  component  parts,  a  sophisticated  algorithm  must  be 
employed.  To  make  this  procedure  as  inexpensive  as  possible, 
two  criteria  must  be  met.  First,  the  number  of  samplings  of 
the  take  off  angle-distance  curve  must  be  minimized  and 
second,  the  tracing  of  each  single  ray,  which  yields  one 
datum  on  the  curve,  must  be  as  fast  as  possible.  The  latter 
problem  can  be  solved  by  using  circular  approximation  ray 
tracing.  The  former  has  also  been  dealt  with  successfully  by 
Hron  (personal  communication)  and  will  not  be  discussed  in 
this  thesis.  One  important  feature  of  Hron' s  process  of 
determining  travel  time  branch  parts  should  be  mentioned. 
When  using  the  linear  method  of  approximating  velocity 
required  by  circular  approximation  ray  tracing,  we  are  in 
essence  introducing  false  interfaces  of  the  second  order. 
Cerveny  and  Pretlova  (1977)  have  examined  this  problem  in 
detail  and  have  drawn  the  conclusion  that  some  degree  of 
smoothing  is  necessary  to  ensure  meaningful  amplitude 
distance  curves.  Hron  has  implemented  a  smoothing  operation 
into  his  program  for  computing  travel  time  branches 
( persona  1  commun i ca  t i on ) . 

One  further  remark  should  be  mentioned  regarding  Figure 
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7.1.  The  infinite  amplitude  predicted  by  asymptotic  ray 
theory  at  the  caustic  C  is  due  to  the  shrinking  of  the  ray 
tube.  It  is  well  Known  that  asymptotic  ray  theory  is  not 
valid  in  the  vicinity  of  points  where  the  phase  function  is 
not  analytic  (see  Cerveny  et  al  1977  for  a  complete 
discussion  of  these  points).  In  such  a  case  it  is  necessary 
to  seek  a  solution  to  the  elastodynamic  equation  in  the  form 
of  another  asymptotic  expansion.  A  method  of  solution  in  the 
vicinity  of  a  caustic,  which  somewhat  resembles  the  ray 
method,  is  developed  in  Appendix  D.  This  alternative  high 
frequency  approximation  will  be  used  in  the  vicinity  of  all 
caustic  points. 


7 . 3  General  mode  1 1 i nq  procedure  for  lateral  1 v  i nhomoqeneous 


medi  a 


A  general  approach  to  producing  ray  theoretical 
synthetic  seismograms  is  given  in  the  flow  chart  of  Figure 
7.2.  The  process  is  for  a  harmonic  source  with  a  single 
predominant  frequency,  but  could  be  adapted  to  sources  with 
more  complicated  frequency  spectra. 

The  efficiency  of  the  method  outlined  in  Figure  7.2 
relies  heavily  upon  the  ray  tracing  procedure  employed. 
Circular  approximat ion. provides  accurate,  yet  rapid  results 
and  is  well  suited  for  implementation. 
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Figure  7.2  General  modelling  procedure  for  laterally 
inhomogeneous  media. 
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Computationally,  the  most  costly  portion  of  the 
modelling  process  is  the  determination  of  individual  travel 
time  branches.  However,  once  computed,  these  results  can  be 
conveniently  stored  on  a  disk  file  in  the  computer  memory 
for  scanning  and  implementation  in  subsequent  parts  of  the 
procedure.  This  feature  also  allows  one  to  determine  the 
effect  on  synthetic  seismograms  of  individual  parts  of 
various  travel  time  branches.  As  well,  a  unique  code  number 
can  be  assigned  to  every  part  of  each  travel  time  branch. 
This  code  number  can  then  be  used  to  identify  events  on  the 
synthetic  traces.  This  is  a  key  advantage  of  ray  methods 
over  full  wave  methods,  such  as  finite  differences,  in  that 
the  user  can  determine  what  portions  of  the  subsurface  have 
been  sampled  in  order  to  produce  a  particular  event. 

The  loop  over  rays  in  Figure  7.2  employs  a  modification 
of  the  numerical  ray  generation  scheme  as  proposed  by  Hron 
(1972).  Since  the  presence  of  lateral  inhomogeneities 
precludes  the  assemblage  of  rays  into  kinematic  and  dynamic 
analogs,  that  portion  of  the  algorithm  has  not  been 
utilized.  Instead  the  rays  are  considered  for  tracing  and 
amplitude  calculations  individually.  The  author  has  compiled 
a  complete  package  of  FORTRAN  programs  employing  the  process 
shown  in  Figure  7.2.  Some  subroutines  were  utilized  which 
were  coded  by  F.  Hron,  namely  the  approximation  of 
interfaces  as  described  in  Appendix  C,  discrimination  of 
travel  time  branches  and  plotting  of  synthetic  seismograms. 
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These  programs  have  been  tested  against  Known  results  for 
simple  media.  Numerical  results  for  two-dimensional 
laterally  inhomogeneous  media  follow  in  the  remaining 
sections  of  this  chapter. 


7 . 4  Cerveny  model 

We  wi  1 1  first  consider  a  vertically  inhomogeneous  model 
utilized  by  Cerveny  and  Zahradnik  (1972)  to  study  the  wave 
field  near  a  caustic.  It  consists  of  an  upper  layer  with 
linear  velocity  gradients  overlying  a  homogeneous  halfspace. 
The  ray  diagram  for  primary  P  waves  and  a  description  of  P 
wave  velocity  in  the  model  are  shown  in  Figure  7.3.  Also 
illustrated  are  the  travel  time  curves  for  these  waves.  The 
threefold  travel  time  branch  consists  of  one  reflected  part 
and  two  refracted  parts  which  form  a  caustic  on  the  surface. 
By  analytical  techniques  ( Brekhovskikh  1960,  section  38,  or 
Choi  1978)  the  location  of  this  caustic  can  be  shown  to  be 
at  120  km  which  is  very  close  to  the  result  obtained  by 
circular  approximation  ray  tracing. 

The  vertical  amplitude  distance  curves  of  the  three 
branch  parts  are  also  given  in  Figure  7.3.  The  method  of 
circular  approximation  was  used  in  all  regions  except  the 
vicinity  of  the  caustic  where  the  Airy  approximation  (see 
Appendix  D)  was  utilized.  The  frequency  of  the  source  was 
taken  as  15  Hz.  The  dotted  line  represents  the  amplitude 
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Figure  7,3  Rays,  travel  times  and  amplitudes  for  the  Cerveny 
model.  The  dotted  curve  was  computed  by  analytical  means 
( Choi  1978) . 
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predicted  by  the  analytical  approach  of  Choi  (1978)  valid  in 
the  neighborhood  of  a  caustic  in  a  vertically  inhomogeneous 
medium.  The  results  match  very  well  over  the  entire  zone  of 
applicability.  Calculation  of  wave  amplitude  by  the  Airy 
function  approach  was  limited  to  those  regions  for  which  the 
difference  in  travel  times  between  rays  lying  on  each  branch 
part  was  less  than  the  time  duration  of  the  source  pulse. 

Figure  7.4  shows  synthetic  seismograms  computed  for  the 
Cerveny  model  in  the  vicinity  of  the  caustic.  The  reduced 
travel  time  branch,  correspondi ng  to  reflected  and  refracted 
paths  has  been  overlaid  onto  the  seismograms.  The  increased 
energy  at  the  caustic  is  clearly  visible  in  these 
sei smograms . 


7 , 5  Bohemi an  Massi f 

The  Bohemian  Massif,  illustrated  and  employed  in 
Chapter  6,  was  used  extensively  by  Psencik  (1972)  to  study 
travel  times  of  primary  P  waves.  This  model  for  the  earth's 
crust  allows  a  full  demonstration  of  the  general  modelling 
process  for  seismic  waves  in  laterally  inhomogeneous  media. 
Figure  7.5  shows  the  ray  diagram  of  reflected  and  refracted 
P  rays  from  a  source  at  155.2  km  for  rays  arriving  at 
prescribed  receiver  points.  For  some  receiver  locations 
there  may  be  up  to  six  arrivals,  although  some  may  be  of 
very  low  amplitude.  The  presence  of  internal  caustics 
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REFLECTED  &  REFRACTED  P  RAYS  IN  THE  BOHEMIAN  MASSIF 
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(corresponding  to  points  where  the  ray  tube  vanishes)  can  be 
ascertained  by  examination  of  Figure  7.5. 

For  reflected  and  refracted  P  waves,  a  computer 
program,  employing  circular  approximation  ray  tracing  and  an 
optimized  search  technique  of  the  distance  -  take  off  angle 
function,  distinguished  11  travel  time  branches.  Some  of 
these  branches  had  up  to  four  parts.  Time-distance  and 
amplitude  distance  curves  of  some  of  the  arrivals  are  shown 
in  Figure  7.6.  Only  the  stronger  arrivals  have  had  their 
amplitude  plotted.  The  Airy  function  approach  for  a 
laterally  inhomogeneous  medium  was  employed  to  compute 
amplitudes  near  two  of  the  caustic  points,  showing  that  some 
of  these  amplitudes  may  be  of  non- trivial  magnitude.  The 
predominant  source  frequency  was  taken  as  TO  Hz. 

Synthetic  seismograms  for  the  Bohemian  Massif  have  been 
constructed  and  are  displayed  in  Figure  7.7.  In  addition  to 
primary  P  waves,  labelled  P1P1,  converted  phases,  PI  SI, 
SI  PI,  and  a  multiply  reflected  wave  (PI)4  were  included. 
Other  rays  were  not  admitted  because  of  either  low  amplitude 
relative  to  P1P1  or  else  very  late  arrival  times. 

One  can  observe  several  phenomena  present  in  the 
synthetic  seismograms  displayed  in  Figure  7.7:  a  strong 
interference  of  multiple  arrivals,  abrupt  appearance  and 
disappearance  of  wavelets,  and  a  great  range  of  amplitudes. 
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SYNTHETIC  SEISMOGRAMS  FOR  THE  BOHEMIAN  MASSIF 
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Figure  7.7  Synthetic  seismograms  for  the  Bohemian  Massif. 
Only  reflected  phases  have  been  interpreted  with  travel 
times.  All  other  arrivals  correspond  to  refracted  branches 
of  PI  PI  or  (PI)4. 
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These  factors  all 
synthetic  traces, 
number  of  internal 
in  these  seismograms 


lead  to  the  variegated  nature  of  the 
The  phase  shift  of  m x/z  where  n  is  the 
caustics  (Choi  1978),  is  also  noticeable 


7 . 6  Conclusions 


While  the  previous  chapter  established  the  accuracy  of 
tracing  rays  and  computing  amplitudes  by  the  circular 
approximation,  this  chapter  has  dealt  with  applications  to 
theoretical  and  realistic  models  of  the  earth's  crust.  The 
numerical  results  presented  support  the  value  of  the  general 
modelling  process  described  in  Figure  7.2.  Clearly, 
asymptotic  ray  theory,  if  properly  employed,  will  play  a  Key 
role  in  efficient  approximate  reconstruct i on  of  geological 
structures  from  field  data. 
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APPENDIX  A 

EQUATIONS  FOR  RAY  INTERSECTION  POINTS 


In  this  appendix  we  shall  derive  explicit  equations  for 
the  location  of  intersection  points  between  the  normal  to 
the  central  ray  and  a  nearby  auxiliary  ray. 

Let  the  equation  of  the  central  ray  be 

(x-a) 2  +  (z-b) 2  =  p2  /  a  i  \ 


where 

a  =  vQcote0/k 
b  =  v  g/k 
p  =  vQ/k  sinQ 

vQ  =  velocity  at  the  beginning  of  the  ray  segment, 
k  =  velocity  gradient  in  the  current  layer, 

0q  =  take-off  angle  at  the  beginning  of  the  ray  segment. 
Equation  A.1  is  given  in  a  Cartesian  coordinate  system  with 
its  origin  at  the  point  where  the  current  ray  segment  begins 
( see  Figure  A.1). 

Let  an  auxiliary  ray  be  given  in  these  coordinates  by 

(X-(A+Xg)  )  2  +  (Z-B)  2  =  R2 


(A. 2) 


' 


,  J 

V 


(X  ,Z)  =z(0 , 0)  Interface 
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Figure  A. 1  Geometry  employed  to  find  ray  intersection 
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where 


A  =  vQcot(j)0/k 
B  =  vQ/k  =.  b 
R  =  v Q/k  sinc{)0 

x  =  lateral  separation  between  starting  points  of  the 

2 

central  and  the  auxiliary  ray.  If  the  ray  segment  under 
consideration  is  the  first  one  emanating  from  the 
source,  then  xg  =  0 . 

We  seek  the  values  of  the  intersection,  (x  ,z  )  ,  of 
the  auxiliary  ray  (equation  A. 2)  and  the  line  M  normal  to 
the  central  ray,  namely, 


z 


mx 


mx  .  +  z  . 
l  i 


(A. 3) 


where 


m  =  (b-zi) / (a-xi) 


Substitution  of  equation  A. 3  into  A. 2  and  introduction  of 


yields 


S*N  +  TXN  +  U  =  ° 


the  solutions  of  which  are 


x 


■N 


( -T  ±  (T2-4SU) ^)/2S 


(A. 4) 


. 

. 


. 
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The  ±  sign  is  taken  according  to  which  solution  gives  an 
intersection  nearest  the  interface  in  the  direction  of  ray 
propagation.  By  substitution,  we  have 


ZN  =  mXN  '  mXi  +  Zi 


(A. 5) 


It  is  also  valuable  to  know  certain  additional 
quantities  concerning  the  auxiliary  ray's  intersection  with 
the  line  normal  to  the  central  ray.  The  travel  time  along 
the  auxiliary  ray  from  (x  , o)  to  (x  ,z  )  is 

o  L\J  IN 


1 

k 


cosh 


-1 


k2(ZN  +  (XN~XS)2)2 
2v0(V0+kZN) 


1 


( A  .  6  ) 


The  length,  &N  ,  from  (x  ,  z  )  to  the  central  ray  as 
measured  along  the  normal  is 


(x  -x  .  )  *  + 
N  i 


(A. 7) 


The  angle 


at  (XN,ZN^ 


that  the  auxi 
i  s 


ary  ray  makes  with  the  vertical 


<j> 


2tan 


_  1 


kt 


N 


tancp  q/ 2 


(A. 8) 


b  i 
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The  degeneration  of 
homogeneous  layers  (k  =  0) 
is  obtained  by  taking 
expressions . 


equations  A.4-A.8  for  cases  of 
or  vertical  incidence  (  0q=°^7T.  ) 
the  limiting  values  of  these 


. 
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APPENDIX  B 

THE  WEBER  FUNCTION  APPROACH  FOR  AMPLITUDE  COMPUTATION 


The  integral  representing  the  vertical  displacement  due 
to  a  harmonic  wave  of  angular  frequency  co  which  has 
penetrated  L  layers,  being  reflected  M  times  from  the  L+1-th 
interface,  and  arrives  at  the  free  surface  at  an  epicentral 
distance  r  is  written  as 

RM(q)a(q)  Hq (krq) e  lkB(q)q  dq 

(B.  1  ) 


A  (r)  = 


ike 


where  we  have  used  the  following  notation 

i  =  f-1 

k  =  10/ a., 

r  ...  plane  wave  reflection  coefficient  for  the  L+1 
interface , 

a  ...  product  of  all  other  plane  wave  coefficients 

correspond!' ng  to  the  wave's  traverse  of  the  medium, 
B  . . .  phase  function  of  the  wave, 
hJ  . . .  Hankel  function  of  the  first  kind  and  zeroth  order. 


After 

1965, 

funct 


suitably 
Cerveny 
on  for  hi 


deforming  the  contour  of  integration  (Cerveny 
and  Ravindra  1971)  and  expanding  the  Hankel 
gh  frequencies,  approximations  may  be  applied 


■ 
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to  the  integral  when  integrating  along  a  portion  of  the 
steepest  descent  contour  given  by 


(1-q  )  2  =  (l-3 


:  ) 


+  pe 


-i.Tr/4 


— oo<p<o° 


(  B  .  2  ) 


The  saddle  point,  x,  is  given  by 


,;/dBX'i 

[Aqj q=x 


=  r 


(  B  .  3  ) 


Physically,  x  is  closely  related  to  the  constant  C  in 
Snell's  Law  by  the  relation  x  =  a^C  =sin01  where  0  is  the 
angle  of  incidence  of  a  P  wave  on  the  first  interface.  For 
computational  purposes,  x=sin01  has  been  chosen  as  a 
parameter  for  the  ray  in  this  thesis.  Ray  parameters 
corresponding  to  critical  refraction  along  the  L+1  interface 
can  be  written  as  x£  =  a./vL+1?  i  =  1,...,L,  where  vL+p  is 
the  velocity  of  the  refracted  segment  of  the  head  wave. 


Due  to  the  singularity  of  R  at  q  =  x£  ,  corresponding 
to  the  critical  point  in  terms  of  the  ray  parameters,  x,  the 
reflection  coefficient  R  can  be  conveniently  decomposed  as 

R(xl)  =  E(xl)  -  D(xl)S(xl) 


where  D  and  E  depend  only  on  even  powers  of  S. 


0  difl 
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Performing  these  operations  on  equation  B.1  yields  an 
expression  for  the  displacement  contribution  of  the 
reflected  wave  near  the  critical  point: 


Vr) 


=  e 


-iu)(t  tR)  °(V  'pM(^  -  MEM_1(xl)D(xl) 


■  *2  2\!s  >'gl(B) 

c;.,xl  -  xl  )  \ 


-  1  i  > 


(  B .  5  ) 


where  the  following  notation  has  been  employed 


z  time  of  arrival  of  the  reflected  wave 

0  r  geometrical  spreading  of  the  reflected  wave, 

1  for  P  type  head  waves 

C  = 


L+l 

LL+1 


for  S  type  head  waves 


x .  =  a1  x/a^, 


l — 1/  •  •  •  /  -i— 1 


6  = 


k (r-r)  \ 
2x3 


2  ^ 


^3  2  ^ 

f  (i-x*;) 2  -  (i-xz)  - 


g1  (B)  =  it 


-iir/4  -p^  j 
:  3-pe  e  dp 


NPi'  NSi  number  of  P  and  S  ray  segments  in  the  i-th  layer, 


h^  thickness  of  the  i-th  layer 


r  =  x 


L  N  . h .  (a^/a?  -  1)  L  N  h,  (a^/b?  -  1) 
3  r  Pi  ill  ,  V  ^  1  1 _ 

l  on  9-3/0  L 


2,2  2,  3/2  ‘  .  ,  2  .,2  2,  3/2 

i=l  (a^/a^  -  x  )  /  i=l  (a1/bi  -  x  ) 


For  regions  beyond  the  critical  distance, 


' 
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r  ?  the  integral  in  equation  B.1  has  contributions  from  a 
circumvention  of  the  branch  cut  given  parametrically  by 


(l-qZr2  =  (1 


2  k  -in/ 4 

x*  )  2  +  pe 


0<p<‘ 


(  B  .  6  ) 


Physically,  this  may  be  interpreted  as  the  head  wave 
contribution  to  the  displacement  on  the  surface. 
Mathematically,  the  omission  of  the  integration  along  this 
path  would  prohibit  a  transformation  of  integration  contours 
since  we  would  end  up  on  a  Riemann  sheet  different  from  that 
at  the  beginning  of  the  contour.  Evaluating  equation  B.1 
along  contour  B.6  and  remembering  that  the  maximum  value  of 
the  integrand  occurs  near  p=0,  gives  a  head  wave 
contribution 


A„(r)  =  Me 

il 


-ico  (t-tR) 


a  (x* )  T  E 

Li 


M-l 


(x* )  c 


X *  (1-x*2) 
8kr  2 (r-r *) 2 


g2  <Y*) 


where 


t 


H 


head  wave  arrival  time, 

head  wave  coef  f  icient  ,=  D(x£)x£ 


(  B  .  7  ) 


f*  =  f  with  x*  substituted  for  x 


y* 


(l-x*2)/ 

- - - - 

(r-r*)  j 


(r-r*) 


g2(y*) 


7iTT/8 

e 


f  h  .  2  „  iu/4  . ,  - 

p^exp(-p  -2e  py*)dp 

0 


■ 


. 
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Near  the  critical  point  we  must  consider  the  summed 
total  response  due  to  the  head  wave  interfering  with  the 
reflected  wave.  Adding  equations  B.5  and  B.7  and  making 
approximations  which  are  legitimate  within  the  interference 
zone,  gives  the  displacement  of  the  interference 
ref lected-head  wave 


A  ( r ) 


-ioo  ( t-tR) 
e 


a(xL) 

Q 


RM (x  ) -MFc 


*2.  % 
( 1 — x  *  1 ) 

x*k (r-r* ) 


,M-1 


(x*) 


i23/4(g1(B)-3! 


i)_2-3/4e-iy 


g0 (y*) 


( B .  8 ) 


where  y  =  (w(t  -t  ))'s  •  A  new  function,  G,  is  introduced  as 

R  H 


G  =  i23/4(g1(B)-e!s)  -  3/4  e'iy  g2(Y*) 


(  B  .  9  ) 


Numerical  investigation  has  shown  that  for  realistic 
sei smologi ca 1  models,  the  three  parameters  8,  y  and  y* 
are  nearly  equal  within  the  interference  zone.  They  can  be 

assigned  a  common  definition  as 

3  for  r<r* 


y 


(0)  (tR-tH)  ) 


for  r>r* 


With  the  introduction  of 
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Y  =  p-ye1Tr//4  in  9-j_(y)  ,  G(y)  can  be  recast  as 

oo 

G(y)  =  -23^4  tt  2  exp  ( 7i7r/8-iy2) 


Y 


exp(-Y2-2yYe17T/4)dY  -  i23/4y^ 


(B. 10) 


This  expression  for  G  can  be  rewritten  with  the  help  of  the 
integral  form  of  the  Weber  function  D  (Magnus  and 


Oberhet t i nger  1949), 


.  .2  .  \ 

Di,  (y  ( i-1 )  )  =7r  ^  2  4  exp  -  ^ - 


H  2  in/M. 

Y  exp(-Y  2yYe  )  dy 


•  oo 


Near  the  critical  point  S-0  in  equation 
formula  B.8  may  be  simplified  further 
eM_1(x*)  =  RM  1(x£)  •  With  these  changes  and 
pretabulated  Weber  function  file,  equation  B.8 
to  rapid  evaluation  on  a  digital  computer. 


B.4  and  thus 
by  setting 
the  use  of  a 
lends  itself 


. 


' 


. 
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APPENDIX  C 

APPROXIMATION  OF  INTERFACES  BY  PIECEWISE  CONTINUOUS 

PARABOLAE 


A  method  of  piecewise  approximation  of  curvilinear 
boundaries  by  a  group  of  smoothly  linked  parabolae  was  first 
suggested  by  Marci novskaya  and  Krasavin  (1968).  Due  to  the 
inaccessibility  of  this  paper,  the  method  is  presented  here 
in  the  modified  form  devised  by  Hron  (personal 
communication)  and  used  in  his  programs  since  1973.  . 


Suppose  that  a  set  of  M  discrete  points,  specified  by 

their  two-dimensional  Cartesian  coordinates  N  =  [L,5m ]  , 

m  mm 

m=1,...,M,  are  given  as  points  lying  on  a  curvilinear 
boundary.  We  may  approximate  this  boundary  with  the  help  of 
I  parabolae 


z (x)  =  (aix+bi) x  +  ci , 


i=l 


I 


(C.  1  ) 


each  of  them  being  defined  over  the  interval 

X£ [Xi,Xi+1l 

As  continuity  conditions  we  impose 

zi(Xi+l)  =  2i+l(XiH! 


dz 


dx 


dz 


i+1 


xi+i 


dx 


xi+i 


21.  .'HIT  51 


'  < 


,  *  ^ 
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The  continuity  of  the  curve  that  passes  through  all  M 
points  could  be  realized  by  solving  the  following  linear 
system  for  the  (i+1)-th  parabola  given  the  i-th  one: 


o 


xi+l 

Xi+1 

1 

ai+l 

zi(xi+i> 

"m 

«a 

1 

bi+l 

Sm 

(C.2) 

-2xi+i 

1 

0 

-  Ci+1 — 

-2aiXi+l  +  biXi+l 

where 

are 

coordinates 

of  the  closest  input 

point 

for  which  Xi+2  =  ^m>Xi+l  •  However,  if  this  system  were 
solved  for  all  M-1  pairs  of  adjacent  points,  considerable 
CPU  time  would  be  consumed  and  the  surface  would  have  a 
corrugated  nature  for  some  types  of  input  data.  The  problem 
of  time  reduction  can  be  vanquished  by  allowing  the  curve  to 
fall  short  of  some  input  points  by  an  arbitrarily  specified 
distance.  For  this  purpose,  equation  C.2  may  be  solved  for 
some  more  distant  point  (iteratively,  beginning  with  point 
M)  and  the  i+1-th  parabola  will  be  accepted  if  for  any  e 


(ai+l^j+bi+l) +  ci+l 

Max ( C j )  -  Min ( 5 . ) 


< 


for  , 


Xi+l<5j<5m 


This  method  is  defeated  if  the  C.  are  scattered  over 


a  large  interval  compared  to  the 


* 


0 
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interval.  This  will  bring  about  the  undesirable 
corrugated  surface.  In  order  to  avoid  this,  we  must  modify 
the  approach.  First,  a  parabola  is  found  that  passes  through 
N£'N£+1'N£+2  by  solving  the  usual  set  of  three  linear 
equations  for  the  coefficients  a^+1 'b£+1 'cz+i  •  Then  the 
two  predetermi ned  parabolae  (a  ,  +  1 ,b ?±1 /C^±1)  will  be 

united  by  an  inserted  parabola  (a^b^c^)  .  The  right 

bound  of  the  interval  of  definition  of  the  1-th  parabola  is 

The 

quantity  Ax  is  selected  iteratively  to  cause  the  relation 


termed  x£+1  and  its  left  bound  is  =  £p-Ax 


XZ<^  Z<XZ+1<^  £+1 


to  hold.  The  coefficients 


a  = 


~2aZ-lXl  +  2a£+lX£+l  b£-l  +  b  i+1 

2(X£+l"X£} 


hZ  2(a£+l  a£)X£+l  +  b£+l 


CZ  (a£  a£+l)X£+l  +  C£+l 


( C  .  3  ) 


and  the  right  bound 


x 


£+1 


2(cm~  Vi*  +  (b«-+i  ~  b<t-i)xa 

2(ai-fam)x!i  +  (b£-i-b£+i) 


( C .  4 ) 


are  the  solutions  of  the  following  systems  of  algebraic 
equations  (with 
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as  a  parameter) 


a&-lX£  +  b£-lX£  +  C£-l  a£X£  +  b£X£  +  C£ 


2a£-lX£  +  b£-l  2a£X£  +  b£ 


a£X£+l  +  b£X£+l  +  C£  a£+lX£+l  +  b£+lX£+l  +  C£+ 1 


2a£X£+l  +  b£  2a£+lX£+l  +  °£+l 
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APPENDIX  D 


RAY  AMPLITUDE  IN  THE  VICINITY  OF  A  CAUSTIC 


Asymptotic  ray  theory,  despite  its  wide  degree  of 
relevance,  is  sever ly  limited  in  the  neighborhood  of  any 
point  where  the  rays  cross  each  other.  Such  a  point  is 
termed  a  caustic.  In  this  appendix  the  el astodynamic 
equation  of  motion  will  be  solved  for  an  arbitrary  ray  in  a 
laterally  inhomogeneous  medium.  The  works  of  previous 
researchers  (Bremmer  1949,  Brekhovskikh  1960,  Ludwig  1966, 
Sato  1969,  Zahradnik  1970,  Chapman  1976,  1978,  1979,  Choi 
1978,  for  example)  have  focused  attention  on  generally 
homogeneous  or  vertically  imhomogeneous  media.  Gazaryan 
(1960)  has  presented  the  solution  to  the  wave  equation  in 
the  vicinity  of  a  caustic  by  a  modification  of  the  ray 
method.  Yanovskaya  (1964)  considered  the  problem  of  a 
solution  to  the  e 1 astodynami c  equation  in  the  vicinity  of  a 
caustic.  Since  these  last  two  papers  are  inaccessible  to 
most  readers,  a  brief  derivation  of  these  results,  with 
modifications  suitable  to  the  geometry  used  for  this  thesis, 
will  be  presented . 

Consider  the  el astodynamic  equation  of  motion  for 
displacement 


L(u) 


(D.  1  ) 


' 


.  •  • 


. 
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where 


L  (u) 


=  (A+u ) V ( V  *u)  + 


iiV23  + 


(V*u)VA  + 


VlixVxu  +  2  (VP  •  V)  u 


in  Cartesian  coordinates  x;i_-x'  x2~^ '  x3_z*  Let  us  impose 

-ia)t 

a  time  dependence  e  .  Following  the  basic  geometry  of 

the  ray  method,  ray  coordinates  (  qi'g2'x  )  are 
introduced  as 

^  . . .  take-off  angle  at  source  measured  from  vertical , 
q2 . . •  azimuthal  take  off  angle  at  source, 

T  ...  phase  function  along  ray. 

We  seek  a  solution  of  D.1  in  the  form 


t)  e 


-iaj  ( t-T ) 


(D.2) 


Asymptotic  ray  theory  states  that  the  following  are 
solutions  to  equation  D.1  (  to  zeroth  order) 


u 


OP 


0 (qirq2)eT 
(Hi^paf1 


( for  P  waves ) 


uos 


ijj  (q  ,q2) 

(H1H2pb)" 


e  cos?  +  e  sin^ 

qi  q2 


( for  S  waves ) 


where 


■ 
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< PrV ...  solve  the  el astodynamic  equation  and  boundary 
conditions  at  each  interface, 
p  .  .  .  volume  densi ty 

a'b  ...  velocities  of  propagation  of  P  and  S  waves 
H1'H2  •  •  •  scale  factors  (Riemann  metrics) 


H 


1 


unit  vectors  at  any  point  along  the  ray, 
angle  between  the  plane  of  incidence  and  the 
direction  £ 


At  a  caustic  H,H2=0  .  We  assume  that  ^=0,  H?/0 

at  this  point.  Physically,  this  excludes  rays  which  focus  at 
a  point.  Furthermore,  we  presume  that  the  rays  propagate  in 
the  (  .q1f t  )  plane  implying 


u(qi,T)  =  u1eg^  +  u3§x 


Substituting  equation  D.2  into  D.1  yields 


N  (u)  ( ico )  2 


(  D  .  3  ) 


where 

N  ( u)  =  (A+y)  (u*Vt)Vt  +  Uu(Vt)2  -  pu 

M  (u)  =  ( A+y  )  (V  t  (V  •  u)  +  V  (u  •  V  x )  )  + 

AA(u*7t)  +  (Vy*u)Vx  +  (Vy»Vx)u 

+  y  (uV2t  +  4^  (Vt)  2) 

d  T 
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We  shall  only  consider  here  the  case  when  the  last  ray 
segment  is  propagating  in  the  compress i ona 1  (or  P)  mode. 
Similar  analysis  for  an  S  wave  is  possible.  For  a  P  wave, 
the  eikonal  equation  is  (Vx)  =  a  .In  what  follows, 
components  along  q1  will  be  indexed  1,  and  along  t 
labelled  3.  Thus  N3  (u)  =  o  ,  ^(u)  =  -pCL-y2)^,  y=b/a 

and  equation  D.3  has  the  form 

L~  (u)  +  iajM_  (u)  =  0 

3  (  D  .  4 ) 


L^(u)  +  ioiM^(u)  -  ioop  ( 1-y  ^ =  0 


(  D  .  5  ) 


Let  us  introduce  a  quantity  n  -  H1H2pa  .  Then  using  the 
definition  of  vector  operators  in  curvilinear  coordinates 


(Spiegel  1959),  we  have 
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[  1  3li  I  H1  3  U1  +  a  3  U3 


a  3t  a  3t  H. 


hi  3qi 


(  D  .  6 ) 
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(  D  .  7  ) 
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(  D  .  9  ) 
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Set  t  =  ?(q  )  as  the  equation  of  the  caustic  surface.  It  is 
also  assumed  that  on  the  caustic  H2^0'  3H1/3t^0  .  Then 


n (q1,x) 

and  therefore 


in 

3t 

T=T 


(T  -  T(q1)  )  +  . .  . 


dq± 


in  3f 

3t  ~  3q.. 

t=t  ^1 


(D. 10) 


We  seek  a  solution  to  equations  D.4  and  D.5  in  the  form 

ux  =  v ( qx )  $1(n) 

u3  =  v(q;L)  $3(n)  (D.  1 1 ) 

Ecjations  D.11  can  be  substituted  into  equations  D.6  to  D.9 
and  subsequently  into  D.4  and  D.5.  Terms  like  VA,Vu,Va 
are  of  a  highly  insignificant  order  of  n  and  can  be 
neglected.  On  the  caustic  n=0  and  the  chief  terms  of 

may  be  determined  only  if  we  retain  the  dominant 

2 

terms.  Clearly,  and  n  $,r  have  the  same  order  in 

n  .  Introducing 
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(D.  12) 
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and  three  operators  of  the  same  order  in  n 

P  ( $)  =  $  +  2ri$  ' 

Q ($)  =  n$ '  -  n2$ '  ' 

R($)  =  $  -  n$'  -  n2$'  ' 


we  have  for  equations  D.4  and  D.5 
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It  is  asssumed  that  co>>i  .  Thus,  in  equation 

2 

disregard  y  P(^1)  since  for  small  n  it  wi 

2  3 

compared  with  MQ ( .$  ) /iooy  n  and  for  large  n 
small  compared  with 


(D. 13) 
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(D. 15) 
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iwn  (1-y2)  yon/3T) 

and  substituted 
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f  3n  YU 

J  \  N  J 


Equation  D. 14  will  be  solved  for 
into  D.15.  In  the  resulting 
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finally  (i-Y2)n<l>{ 


and 


is  neglected  for  the  same  reasons  as  y  ?($-]_)  •  After  these 

processes  we  are  left  with 
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(D. 16) 


In  equation  D.16  the  R($3)  term  is  neglected  with 
respect  to  Q($n)  .  Thus,  the  system  D.14-D.15  reduces  to 

-  7^3  q($3)+  !t'2^r(V  - 

103Y  _ 

+  P($  )  -  N(l-Y2)^1'/On/3T)Z  = 


0 


(D. 17) 
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9 

<- 


m  n ( l— )  io)n(l-y  )$, 

— rr  °<V  +  —  f  2  «;  +  1 


io)y  ri 


( 3ri/3  t  ) 


2  ‘3  ( 3  n/3  t ) 
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(D. 18) 


These  last  two  equations  can  be  expressed  equivalently  as 
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N ( 1-y  ; 


io)  r) 


y  Q($3)  +  P  ( ^ 3 )  - 
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(3n/3  t) 


=  0 


(D. 19) 


n  (l-y2)  $1  io)nd-y  H 


( 3  n/3  t ) 
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(  3 ri/3  t  ) 


=  0 


( D . 20  ) 


From  D.20  we  find  that 


Q($3) 
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(  D  .  2 1  ) 


Substitution  into  D.19  yields 
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io)y  ri 


2“3  Q(^3)  +  P($3)  =  0 


( D . 22 ) 


For  rays  approaching  the  caustic  M<0,  while  for  rays 
leaving  the  caustic  M>0.  These  two  cases  must  be  considered 
separately  when  solving  equation  D.22.  We  define 


. 


, 


163 


,3  I 
£  =  a  I  n 


3 


( D . 23 ) 


2  ,  , 

where  a  =  u)y  /3|M|  .  Then  equation  D .  22  may  be  recast  as 

$2  (1  +  6iari^)  ±  3ar)3$2  =  0 


where  the  upper  (lower)  sign  is  for  the  ray  approaching 
(leaving)  the  caustic.  Substituting 


yields  an  equation  which  is  a  solution  in  the  form  of  the 
Hankel  functions.  By  increasing  the  frequency  to  we  can 
make  the  argument  of  the  Hankel  function  large  in  spite  of 
n  being  small.  If  one  employs  the  asymptotic  expansion  for 
the  Hankel  functions  (Abramowitz  and  Stegun  1965),  the 
zeroth  order  asymptotic  ray  theory  result  must  be  found, 
viz. =  n  2  .  This  allows  us  to  choose  suitable 
normalizing  factors.  Such  a  process  as  this  also  enables  us 
to  define  V(q^)  in  order  to  make  equation  D.11  consistent 
with  the  zeroth  order  theory.  Finally,  we  have 


( D . 24 ) 
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( D . 25 ) 
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V(qj_) 


{vopo) 
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K 
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v  p  cos9  \ 

\  v'  p  cosG  ^ 


R  . 
D 


( D . 26 ) 


where  the  superscripts  "ap"  and  "1e"  denote  rays  approching 
and  leaving  the  caustic.  vq'po  are  wave  velocity  and 
density  at  the  source,  K  is  the  total  number  of  ray 
segments,  Rj  is  the  appropriate  plane  wave  interface 
coefficient  which  ends  the  j-th  ray  segment  v,p,0  are 
the  wave  velocity,  density  and  acute  angle  with  the 
interface  normal  on  the  non-incident  side  of  any  interface. 
The  corresponding  properties  on  the  incident  side  of  the 
interface  are  given  by  v,p,0 


The  funct ions 


may  oe  found 


by  equation 


D.20.  They  are 
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( D  .  27 ) 
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( D . 28 ) 


Due  to  the  presence  of  the  frequency  factor, 
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is  smaller  than  $3  .  This  is  related  to  the  fact 

that  describes  the  transverse  component  of  the  P  wave, 

which  is  zero  in  zeroth  ray  theory. 

The  equations  developed  so  far  have  been  written  in  the 
curvilinear  ray  coordinates.  The  computations  are  performed 
in  Cartesian  coordinates  and  to  this  end  the  appropriate 
conversions  must  be  made. 

Consider  the  point  0  where  the  caustic  surface 
intersects  the  surface  on  which  the  receiver  lies.  Place  at 
0  the  origin  of  a  Cartesian  coordinate  system  (x,y,z)  where 
x  points  normal  to  the  caustic  and  z  points  tangent  to  the 
caustic  and  into  the  earth.  Let  M  be  a  point  at  a  distance 
v  from  0  measured  along  the  x  axis  and  v>0  on  the  lit 
side  of  the  caustic.  Gazaryan  (1960)  has  demonstrated  that 


q  (M) 


(D.29) 


T  (M)  =  T  (  0)  + 


( D . 30 ) 


where  the  upper  (lower)  sign  is  for  the  ray  approaching 
(leaving)  the  caustic.  Since 
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n(M)  =  n  ( 0 )  +  (  (a-|_(M)  -■  ot1  ( 0 ))  +  ... 


ll 
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and  n  ( o ) =  o  ,  it  follows  that 
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(D.31 ) 


We  must  express  the  unit  vectors  ex,ez 
a  •  Expand  e_(M)  as 
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Taking  D.29  into  account  and  also  that  e^(0)  =  -e^  we  have 
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Analogously, 


( D . 32 ) 
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( D . 33 ) 


The  upper  (lower)  sign  in  these  equations  corresponds  to  the 
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ray  approaching  and  leaving  the  caustic.  Components  of 
displacement  in  the  x  and  z  coordinates  can  be  written  as 
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( D . 34 ) 
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( D . 35 ) 


In  equation  D.35  the  term  involving  ^  decreases  rapidly 
with  increasing  v  compared  to  the  $3  term  and  can  thus 
be  neglected.  Expand  the  variables  n  and  t  in  equations 
D.34  and  D.35  and  make  use  of  the  definition  of  the  Airy 
function  in  terms  of  Hankel  functions  to  arrive  at 
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( D . 36 ) 


3ttV1/6 

Uz  =  2v(c^i)  exp  (-iw  (t-f ) +3i7T/4)  Ai  (-Coo2//3’.;)  (D.37) 


where 
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3n  N 1/3 


3  T 
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In  order  that  these  equations  can  be  employed  in  the 
ray  programs  proposed  by  this  thesis,  x  and  z  must  be 
oriented  in  the  standard  direction,  viz.,  z  pointing 
vertically  downward  and  x  pointing  to  the  right  of  our 
profile.  This  is  a  simple  conversion  since  the  angle  0  , 
the  caustic  makes  with  the  standard  z  direction  is  found  by 
tracing  a  ray  which  ends  up  on  the  receiver  surface  at  the 
point  of  intersection  between  that  surface  and  the  caustic. 
For  computational  purposes  we  shall  neglect  the  transverse 

component  u  due  to  the  dependence  as  compared  to 

1  /6 

the  compressiona  1  component  which  varies  as 

Finally,  after  further  simplifications,  the 
displacement  for  P  waves  arriving  at  a  receiver  in  the 
vicinity  of  a  caustic  can  be  written  as 


( D . 38 ) 
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x*slocation  where  caustic  intersects  receiver  surface, 
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Similarly,  for  S  waves  arriving  at  the  receiver,  we 
obtain 
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( D . 39  ) 


where 


D  = 


"  9  H  \  2 


,4  9t 
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92x 
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a'b/p  =■  P  wave  velocity,  S  wave  velocity  and  volume  density 

—y 

at  the  point  x* 

All  quantities  in  equations  D.38  and  D.39  are  easy  to 
evaluate  by  the  circular  approximation  method  of  ray 
tracing,  which  renders  the  expressions  very  useful  in 
computing  amplitude-distance  curves  or  synthetic 
seismograms.  As  pointed  out  by  Choi  (1978)  or  Choy  and 
Richards  (1975)  a  tt/2  phase  shift  is  attributed  to  the 
seismic  pulse  each  time  the  ray  passes  through  an  internal 
caustic.  Circular  approximation  ray  tracing  as  developed  in 
Chapter  5  is  well  suited  to  determining  the  location  and 
number  of  these  caustics.  By  simultaneously  tracing  two  rays 
which  are  separated  by  a  small  amount  in  takeoff  angle,  one 
locates  all  intersection  points  of  the  various  circular  arcs 
comprising  the  ray  paths.  Within  each  triangular  region  of 
the  model  that  the  rays  traverse,  analytical  expressions  can 
be  written  for  the  intersection  of  the  two  circular 
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segments.  If  such  a  point  or  points  is  contained  within  the 
bounds  of  the  triangular  region,  then  the  rays  have  passed 
through  an  internal  caustic.  Otherwise,  this  side-by-side 
ray  tracing  process  continues  into  the  next  triangular 
region . 
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